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Introduction 

This  is  the  final  report,  covering  the  period  May  1,  2005  through  August  31,  2008.  In  this 
report,  we  discuss  our  progress  on  nonlinear  systems  analysis,  based  on  Lyapunov  and  dissi¬ 
pation  inequalities,  using  sum-of-squares  (SOS)  decompositions  to  verify  set  containments. 

We  have  improved  our  ability  to  analyze  uncertain  system  dynamics,  including  nonaffine 
parametric  uncertainty  as  well  as  unmodeled  dynamics.  This  entailed  polytopic  covering 
methods  for  graphs  of  vector-valued  polynomial  functions,  and  local  small-gain  theorems. 
We  use  simulation  as  a  key  step  in  aiding  the  nonconvex  search  for  proofs  (Lyapunov  func¬ 
tions)  and  proof  certificates  (multipliers).  Some  aspects  of  the  calculation  are  trivially  par- 
allelizable,  and  we  have  employed  a  9- machine  cluster  to  speed-up  the  analysis  of  uncertain 
systems.  We  also  began  more  detailed  study  of  systems  with  marginally  stable  linearizations 
(ie.,  adaptive  systems).  Finally,  we  made  precise  our  claim  that  these  techniques  represent 
a  quantitative  and  definitive  improvement  over  linearized  analysis. 


Notation:  R[x]  represents  the  set  of  polynomials  in  x  with  real  coefficients.  For  7r  €  R[xj, 

d(n)  denotes  the  degree  of  n.  The  subset  E[x]  :=  {tt\  + - 1-  7r^  :  7Ti ,  •  •  •  ,  7rm  €  R[x]}  is 

the  set  of  SOS  polynomials. 

Uncertain  Systems  Analysis 

Uncertainty  in  the  vector  field  includes  non-affine  dependence  on  a  parameter  vector  8  lying 
in  a  polytope  A,  namely 
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x(t)  =  fo(x(t))  +  ^<5i/i(x(t))  +  ^2gj(8)fm+j(x(t)), 
i=  1  7=1 
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and  unmodeled  dynamics, 


x(t)  =  + 

z  =  h(x) 
w  =  $(2) 

Here  $  represents  unmodeled  dynamics,  for  example,  finite-dimensional,  linear,  time-invariant 
operators  with  specified  upper  bound  on  induced  £2  norm,  for  example  ||$||£2_£  <  1. 

The  tools  we  have  developed  to  address  this  problem  are 

1.  region-of-attraction  analysis  for  systems  with  affine  parameter  uncertainty  using  a  sin¬ 
gle  Lyapunov  function 

2.  local  induced  £2  —*  £2  gain  analysis  for  systems  with  affine  parameter  uncertainty, 
using  a  single  polynomial  storage  function 

3.  covering,  with  a  polytope,  the  graph  of  vector- valued  polynomial  function  over  a  poly- 
topic  domain, 

4.  informal  branch- and-bound 

5.  local  small  gain  theorems. 

Next,  we  illustrate  the  calculations  that  are  possible  with  these  methods. 

Controlled  short  period  aircraft  dynamics,  parametric  uncertainty 

We  apply  the  robust  ROA  analysis  for  uncertain  controlled  short  period  aircraft  dynamics 
(see  the  project  website  for  the  parameters  used  in  the  model) 


’  Coi(£p)  +  <5iCn(£p)  +  £i<73l(*p)  ' 

(Jxp  +  6n  +  6i2<5i 

Xp  = 

<?02(Zp)  +  Si  £[2£p  +  <52<722(Zp) 

+ 

&2i  +  622^2 

0 

where  xp  =  [xi  £2  £3]T  ,  £1,  £21  and  £3  denote  the  pitch  rate,  the  angle  of  attack,  and  the 
pitch  angle,  respectively,  coi  and  cn  are  cubic  polynomials.  <702,  <722,  and  <731  are  quadratic 
polynomials,  £12  and  4  are  vectors  in  7 Z3,  bn,  612, 621  >  and  622  G  and  u,  the  elevator  deflec¬ 
tion,  is  the  control  input.  Variations  in  the  center  of  gravity  in  the  longitudinal  direction  are 
modeled  by  5\  6  [0.99,2.05]  and  variations  in  the  mass  are  modeled  62  6  [—0.1, 0.1].  Note 
that  the  parametric  uncertainty  includes  one  nonaffine  term  (ie.,  ^).  The  control  input  is 
determined  by  £4  =  — 0.864yi  H — 0.321y2  and  u  =  2£4,  where  £4  is  the  controller  state  and 
the  plant  output  y  =  [£1  £3]r .  Define  x  :=  [: rj  £4]  and  the  shape  factor  p(x)  :=  xTx. 
We  applied  a  branch-and-bound  type  procedure  with  d(V)  =  2  and  d(V)  =  4  on  a  9- 
processor  computer  cluster:  after  the  first  B&.B  iteration,  the  cell  with  the  smallest  lower 
bound  is  subdivided  into  3  subcells  and  cells  with  2-nd,  3-rd,  and  4-th  smallest  lower  bounds 
are  sub-divided  into  2  subcells.  Fig.  1  shows  the  lower  bounds  and  upper  bounds.  Note 
that  quadratic  Lyapunov  functions  (several,  as  different  Lyapunov  functions  are  employed 
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Figure  1:  Lower  bounds  for  ^  with  d(V)  =  2  (solid  black  with  “x”)  and  d(V)  =  4  (solid  blue 
curve  with  “o”)  and  0™  (solid  red  with  “o”)  computed  at  the  centers  of  the  cells  generated  by  the 
B&.B  Algorithm  for  the  d(V )  =  4  run.  Dashed  curves  are  for  (computed  values  of)  0^}  where 
6  is  the  center  of  the  cell  with  the  smallest  lower  bound  at  the  corresponding  step  of  the  B&B 
refinement  procedure  for  d(V)  =  2  (dashed  black  with  “x”)  and  d{V)  =  4  (dashed  blue  with  “o"). 

in  different  cells  across  the  parameter  space)  certify  that  all  initial  conditions  Xq  e  R4  satis¬ 
fying  XqXo  <  5.4  are  in  the  region-of-attraction.  Likewise,  a  collection  of  quaxtic  Lyapunov 
functions  certify  that  all  initial  conditions  x0  €  R4  satisfying  XqX0  <  7.8  are  in  the  region- 
of-attraction.  The  smallest  value  of  p  attained  on  divergent  trajectories,  [3nc,  is  8.6  and 
obtained  for  (<5j,  62)  =  (2.039,  —0.099)  and  the  initial  condition  (0.17, 2.65,  —0.10, 1.24). 


Aircraft  dynamics,  parametric  uncertainty  and  unmodeled  dynamics 

Next,  consider  the  same  system  with  additional  unmodeled  dynamics  at  the  plant  input,  as 
shown  in  Figure  2. 


Figure  2:  Controlled  short  period  aircraft  dynamics  with  unmodeled  dynamics  (5P  := 
(<5i,  <$2))- 

The  assumption  is  that  $  is  any  stable,  linear  time-invariant  (this  could  be  relaxed)  operator, 
with  induced  £2  norm  less  than  1.  We  again  repeat  the  analysis  using  both  quadratic  and 
quartic  storage  functions,  which  locally  certify  bounds  on  the  gain  (with  $  removed)  from 
w  to  z  (recall  that  this  system  is  not  globally  stable)  in  the  presence  of  the  parametric 
uncertainty.  These  local  £2  gains  are  used  in  conjunction  with  a  local  small-gain  theorem  to 
yield  results  such  as: 

•  (using  quadratic  storage  functions)  For  all  (finite-dimensional,  linear,  time-invariant) 
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$,  satisfying  ||^>||jC2_iC2  <  1,  assuming  that  the  initial  condition  of  $  is  0,  then  all 
plant/controller  initial  conditions  x0  G  R4  satisfying  XqXq  <  2.4,  are  in  the  robust 
region-of-attraction. 

•  (using  quartic  storage  functions)  For  all  (finite-dimensional,  linear,  time-invariant) 
$,  satisfying  ||$||£a_£a  <  1,  assuming  that  the  initial  condition  of  $  is  0,  then  all 
plant /controller  initial  conditions  x0  G  R4  satisfying  XqXq  <  4.1,  are  in  the  robust 
region-of-at  tract  ion . 

In  conclusion,  we  have  established  provable  and  certifiable  inner  estimates  of  the  region-of- 
attraction  of  a  (nominally)  4-state  nonlinear  system  with  both  parametric  and  unmodeled 
dynamics  uncertainty.  The  informal  use  of  branch-and-bound  in  the  parametric  uncertainty 
space  was  handled  efficiently  using  a  small-scale  parallel  cluster  of  9  machines. 


Viewing  our  approach  as  quantitative  extension  of  linearized  analysis 

Practical  nonlinear  analysis  often  couples  extensive  nonlinear  simulation  with  extensive  lin¬ 
earized  analysis  (such  as  stability,  stability  margins.  Bode  plots  of  linearized  I/O  maps,  etc). 
Here  we  show  that  common  linearized  analysis  techniques  can  be  rigorously  quantified  using 
the  SOS  approaches.  The  next  lemma  is  key  to  these  derivations. 

Lemma:  Let  z(z)  be  a  vector  of  all  monomials  of  degree  2  with  no  repetition.  Let  Q  = 
QT  >-  0.  There  exists  a  positive  definite  matrix  H  =  HT  such  that  xTxxTQx  =  z(x)J  Hz{x ). 

SOS  region-of-attraction  analysis:  For  the  autonomous  system  x  =  /(x),  and  a  positive- 
definite  function  p,  if  there  exist  positive-definite,  radially  unbounded  function  Zj,  positive- 
definite  function  Z2,  SOS  polynomials  Si,  s2  and  S3,  a  polynomial  function  V,  and  positive 
constants  7  and  (3  such  that  V(0)  =  0  and 

V-Zi<=E[x],  (la) 

-  ft/3  -  p)«,  +  (V  -  7)]  €  E[*]f  (lb) 

-  [(7  -  V)s2  +  W/s3  +  Z2]  e  E[x],  (lc) 

then  {x  :  p(x)  <  (3}  C  {x  :  V(x)  <  7}  =:  it,  and  for  all  x(0)  G  S2,  the  solution  satisfies 
x(t)  G  i 2  Vf  and  lim^oo  x(Z)  =  0. 

Now,  consider  the  system  x(t)  =  Ax(t)+ /23(x(Z))  where  A  is  Hurwitz,  and  /2 3  is  a  polynomial 
with  quadratic  and  cubic  terms.  For  any  quadratic,  positive-definite  Z 1,  Z2  and  p,  the  equations 
(1)  are  feasible  using  quadratic  V,  constant  si,  S3  and  quadratic  s2  (suboptimal,  feasable 
values  are  easily  determined  from  A  and  /2 3).  Consequently,  if  the  local  stability  of  an 
equilibrium  point  of  a  system  with  a  cubic  vector  field  is  decidable  using  linearized  analysis, 
then  the  SOS  region-of-attraction  analysis  will  always  yield  a  quantitative,  certified,  inner 
estimate  of  the  region  of  attraction. 

SOS  £2  gain  analysis:  For  the  driven  system  x  =  /(x,  w),z  =  h(x),  if  there  exist  positive- 
definite,  radially  unbounded  function  Zi,  SOS  polynomial  sj,  a  polynomial  function  V,  and 
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positive  constants  7  and  R  such  that  V(0)  =  0  and 


v  -  h  e  £[*], 


(R2  —  V)si  +  W  /  —  WTW  +  ^hTh 


€  E[x], 


(2a) 

(2b) 


then  for  all  w ,  with  ||ty||2T  <  R,  the  solution  from  x(0)  =  0  satisfies  V(x(t))  <  R2  and 

WzWi,t  <l\\M\w 

Now  consider  the  system  x  =  Ax  +  /2(x)  +  fs{x)  +  [B  +  <7i(x)]  w  and  2  =  Cx  +  h2(x)  where 
gi  is  purely  linear,  /2  and  h2  are  purely  quadratic,  and  /3  is  purely  cubic. 

Suppose  the  linearization  (A,  B  and  C)  has  A  Hurwitz,  and  ||C(s/  —  >1)_1B||00  <  7.  For 
any  quadratic,  positive-definite  l\,  there  exist  R  >  0  such  that  the  equations  (2)  axe  feasible 
(possibly  after  scaling  the  state  coordinates,  x  <—  ax,  by  a  computable  scalar  a  >  0)  using 
quadratic  V,  and  quadratic  Si  (suboptimal.  feasable  values  are  easily  determined  from  A , 
/2,  fzi  etc.).  Consequently,  if  the  local  input/output  gain  around  an  equilibrium  point  of  a 
system  with  a  cubic  vector  field  is  bounded  using  linearized  analysis,  then  the  SOS  £2  gam 
analysis  will  always  yield  a  quantitative,  certified,  ball  of  disturbances  such  that  the  same 
gain  bound  holds  for  the  nonlinear  system. 
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Robust  Region-of-Attraction  Estimation 

Ufuk  Topcu,  Andrew  Packard,  Peter  Seiler,  and  Gary  Balas 


Abstract 

We  propose  a  method  to  compute  invariant  subsets  of  the  region-of-attraction  for  the  asymptoti¬ 
cally  stable  equilibrium  points  of  polynomial  dynamical  systems  with  bounded  parametric  uncertainty. 
Parameter-independent  Lyapunov  functions  are  used  to  characterize  invariant  subsets  of  the  robust  region- 
of-attraction.  A  branch-and-bound  type  refinement  procedure  is  implemented  to  reduce  the  conservatism. 
We  demonstrate  the  method  on  an  example  from  the  literature  and  uncertain  controlled  short  period 
aircraft  dynamics. 


I.  INTRODUCTION 

We  consider  the  problem  of  computing  invariant  subsets  of  the  region-of-attraction  (ROA)  for 
systems  with  polynomial  vector  fields  and  bounded  parametric  uncertainty.  Since  computing  the 
exact  ROA,  even  for  systems  with  known  dynamics,  is  hard,  research  has  focused  on  determining 
Lyapunov  functions  whose  sublevel  sets  characterize  invariant  subsets  of  the  ROA  [8],  [9],  [19]. 
Recent  advances  in  polynomial  optimization  based  on  sum-of-squares  (SOS)  relaxations  [12]  are 
utilized  to  determine  invariant  subsets  of  the  ROA  for  systems  with  known  polynomial  and/or 
rational  dynamics  solving  optimization  problems  with  matrix  inequality  constraints  [21],  [15], 
[7],  [14],  [17].  The  literature  on  ROA  analysis  for  systems  with  uncertain  dynamics  includes 
a  generalization  of  Zubov’s  method  [4]  and  an  iterative  algorithm  that  asymptotically  gives 
the  robust  ROA  for  systems  with  time-varying  perturbations  [11].  Systems  with  parametric 
uncertainties  are  considered  in  [5],  [13],  [18].  The  focus  in  [5]  is  on  computing  the  largest 
sublevel  set  of  a  given  Lyapunov  function  that  can  be  certified  to  be  an  invariant  subset  of 
the  ROA.  References  [13],  [5]  propose  parameter-dependent  Lyapunov  functions  which  lead  to 
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potentially  less  conservative  estimate  of  the  ROA  compared  to  parameter-independent  Lyapunov 
functions  at  the  expense  of  increased  computational  complexity. 

This  paper  follows  [16],  using  bilinear  sum-of-squares  optimization  to  determine  invariant 
subsets  of  the  robust  ROA.  The  differences  lie  in  the  allowed  uncertain  parameter  dependence 
and  the  class  of  Lyapunov  functions.  The  approach  in  [16]  employs  parameter-independent 
Lyapunov  functions  for  systems  whose  vector  field  depends  affinely  on  uncertain  parameters 
known  to  lie  in  a  given  polytope.  This  is  reminiscent  of  quadratic  stability  analysis  [3],  where  a 
single  quadratic  Lyapunov  function  certifies  the  stability  of  an  entire  family  of  uncertain  linear 
systems,  usually  described  by  a  polytope  of  linear  vector  fields.  Of  course,  in  both  cases,  using 
a  common  Lyapunov  function  tends  to  yield  conservative  results.  Additionally,  the  restriction  to 
polytopes  of  vector  fields  is  undesirable.  This  paper  partially  alleviates  both  of  these  limitations. 
First,  vector  fields  are  allowed  to  depend  affinely  on  polynomial  functions  of  the  uncertain 
parameters,  and  we  develop  methods  to  cover  these  with  a  polytope  of  vector  fields  (so  that 
[16]  applies).  Additionally,  we  propose  a  branch-and-bound  type  refinement  procedure  [10]  to 
partition  the  uncertainty  set  and  compute  a  different  parameter-independent  Lyapunov  function 
for  each  cell,  hence  implicitly  using  piecewise  constant  (across  uncertainty  space),  yet  parameter- 
dependent  Lyapunov  functions.  In  fact,  in  robustness  analysis  involving  time-invariant  unknown 
parameters,  it  is  common,  [2],  [22],  to  combine  easily-computable  sufficient  conditions  with 
branch-and-bound  strategies,  often  yielding  improved  analysis  results. 

An  alternate  for  the  conservativeness  of  parameter-independent  Lyapunov  functions  is  using 
polynomially  parameter-dependent  Lyapunov  functions  as  proposed  in  [5],  [13].  Although  SOS 
optimization  can  be  used  with  parameter-dependent  Lyapunov  functions,  the  ensuing  optimiza¬ 
tion  problem  is  challenging  because  uncertain  parameters  are  treated  as  additional  independent 
variables  in  the  SOS  conditions,  which  can  greatly  affect  the  size  of  the  semidefinite  programs. 
Moreover,  choosing  a  suitable  and  effective  polynomially  parameter-dependent  basis  for  the 
Lyapunov  function  is  not  intuitive. 

Finally,  we  remark  that  the  methodology  based  on  the  branch-and-bound  algorithm,  applied 
to  robust  region-of-attraction  analysis  here,  is  generally  applicable  to  local  reachability  and  gain 
analysis  of  systems  with  parametric  uncertainty. 

Notation:  R[x]  represents  the  set  of  polynomials  in  x  with  real  coefficients.  For  ir  €  R[x], 
d{i r)  denotes  the  degree  of  n.  The  subset  E[x]  :=  {nj  H - b  7r^  :  trj,  •  •  •  ,irm  €  R[x]}  is 


the  set  of  SOS  polynomials.  For  77  e  1Z  and  g  :  lZn  — >  TZ,  the  77-sublevel  set  ilg^  of  g  is 
defined  as  il9tV  :=  {x  €  lZn  :  g(x)  <  77}.  In  several  places,  a  relationship  between  an  algebraic 
condition  on  real  variables  and  state  properties  of  a  dynamical  system  is  claimed,  often  using  the 
same  symbol  for  a  particular  real  variable  in  the  algebraic  statement  as  well  as  the  state  of  the 
dynamical  system.  This  could  be  a  source  of  confusion,  so  care  on  the  reader’s  part  is  required. 

II.  ESTIMATION  OF  THE  ROBUST  ROA  OF  SYSTEMS  WITH  PARAMETRIC  UNCERTAINTY 
Consider  the  system  governed  by 

x(t)  =  f(x(t),  6),  (1) 

where  5  €  A  C  7Zm  is  the  vector  of  unknown  parameters  and  A  is  a  known  bounded  polytope. 
For  each  6  €  A,  /(•,£)  :  lZn  — ►  lZn  is  locally  Lipschitz  and  satisfies  f(0,6)  =  0.  The 
robust  ROA  is  the  intersection  of  the  ROAs  for  all  systems  governed  by  (1),  i.e.,  fl,5<=A  {xo  € 
7Zn  :  limf_00  xo,  6)  =  0},  where  </?(t;xo,  S)  denotes  the  solution  of  (1)  at  time  t  with 
initial  condition  Xo  and  fixed  parameter  value  5  €  A.  Trivial  extensions  of  results  found  in 
classic  textbooks  [20]  show  that  sublevel  sets  of  appropriate  Lyapunov  functions  are  invariant 
subsets  of  the  robust  ROA.  For  a  subset  D  C  7 Zm  and  continuously  differentiable  function  V. 

define  Md,v  ■=  PlieD  {x  '■  W(x)/( x,S)  <  0}. 

Proposition  2.1:  If  there  exist  7  >  0  and  a  continuously  differentiable  V  :  lZn  — >  1Z  such  that 


V(0)  =  0  and  V{x)  >  0  for  all  x  ±  0,  (2) 

is  bounded,  and  (3) 

^V,-y\{0}  C  Ma.,V,  (4) 

hold,  then  for  all  xo  €  Civ,- y  and  for  all  S  €  A,  </?(f;xo,  S)  exists,  satisfies  </?(f:xo,  <5)  €  Qy.-r  for 
ail  t  >  0,  and  lim^oo  *p{t\  xq,  S )  =  0,  i.e.,  ily,- y  is  an  invariant  subset  of  the  robust  ROA.  < 


We  now  restrict  our  attention  to  a  special  case,  where  the  dependence  of  /  on  6  is  affine,  to 
obtain  conditions  equivalent  to  (4)  for  this  special  case  and  suitable  for  numerical  verification  (a 
generalization  to  polynomial  dependence  on  6  is  treated  in  section  III).  Assume  that  the  vector 
field  in  (1)  is  in  the  form 

m 

i(t)  = /«(*«)  + £>/<(*»), 


(5) 


where  /0,  /i  > • • • i fm  :  7ln  — *■  7ln  are  known  locally  Lipschitz  functions  and  satisfy  ffO)  =  0 
for  i  =  0, 1, . . . ,  m,  and  5  €.  A.  Further,  denote  the  set  of  vertices  (extreme  points)  of  A  by  £&. 
Then,  the  following  follows  from  the  affine  dependence  of  the  vector  field  on  5,  [16]. 

Proposition  2.2:  For  the  vector  field  in  (5)  and  and  continuously  differentiable  function  V, 
the  equality  Ma,v  =  holds.  « 

Consequently,  for  any  continuously  differentiable  function  V  satisfying  (2),  (3),  and 

SV,7\{0}  c  Ms^y,  (6) 

the  sublevel  set  Ovvy  is  30  invariant  subset  of  the  robust  ROA.  In  order  to  enlarge  the  computed 
invariant  subset  of  the  robust  ROA  by  choice  of  V ,  we  introduce  a  fixed,  positive  definite,  convex 
function  p,  called  the  analysis  shape  factor  and  maximize  (3  while  imposing  the  constraints  (2)- 
(3),  (6),  and  :=  { x  G  lZn  :  p{x)  <  (3)  C  This  is  written  as  an  optimization  problem, 

0a‘(V)  :=  ^  jTiax  o  f3  subject  to  (2),  (3),  (6),  and  ilPi0  C  f lv„.  (7) 

Here,  V  denotes  the  set  of  candidate  Lyapunov  functions  over  which  the  maximum  is  computed, 
for  example  all  continuously  differentiable  functions.  In  practice,  p  is  problem-dependent,  and 
chosen  by  the  analyst.  Since  the  form  of  the  certified  inner  estimate  of  the  robust  region-of- 
attraction  is  a  sublevel  set  of  p,  the  sublevel  sets  of  p  should  be  well-understood  (for  in  high- 
dimensions  they  cannot  be  visualized),  and  should  reflect  directionality/scaling  information  that 
the  analyst  is  interested  in  learning  with  regard  to  the  robust  region-of-attraction.  In  order  to 
relax  the  problem  in  (7)  to  a  SOS  programming  problem,  we  require  /o,  fi,  •  •  • ,  fm  and  p  to  be 
polynomials  and  restrict  V  to  be  a  polynomial  in  x  of  fixed  degree.  Further,  we  use  generalizations 
of  the  S-procedure  [16]  to  obtain  sufficient  conditions  for  the  set  containment  constraints  in  (7) 
and  the  well-known  SOS  sufficient  condition  for  polynomial  nonnegativity  [12]:  if  n  e  E[x], 
then  7 r  is  nonnegative. 

Let  Vpoiy  C  V,  Si,  S2,  and  S3  be  prescribed  finite-dimensional  subsets  of  R[ar],  and  denote 
S  =  (Si,S2, S3).  For  a  polytopic  subset  D  of  A  and  positive  definite  polynomials  and  /2 


(typically  U{x)  =  €iXTx  with  small  scalars  e*),  define  /3D(Vpoiy,S)  as 


si  e  E[x],  s2<5  €  E[xj,  s3<5  e  E[x],  for  all  5  6  SD , 

0>o,  7>0,  1^(0)  =  0,  V"  6  Vpoly,  V  -lxe  E[x], 

~  [(P  -  P)s  1  +  (V  —  7)]  €  E[x],  and 
-  [(7  -  V)S26  +  W (/o  +  $ifi)s36  +  k)  €  E[x],  for  all  8  €  SD. 


(8d) 


(8b) 


(8a) 


(8c) 


The  feasibility  of  the  constraints  in  (8)  is  sufficient  for  the  feasibility  of  the  constraints  in  (7). 


Therefore,  /?aOW$)  <  ^(V). 


The  optimization  in  (8)  is  naturally  converted  to  a  bilinear  semidefinite  program  (SDP),  with  3 
“types”  of  decision  variables:  the  free  parameters  in  V,  the  free  parameters  in  the  s  polynomials, 
and  the  free  parameters  introduced  by  the  SOS  constraints.  The  SDP  is  bilinear  in  the  free 
parameters  in  V  and  multipliers  s,  as  evidenced  by  the  product  terms  (e.g.  Vs2,5,  VV/s3,5,  etc). 
We  have  made  significant  pragmatic  progress  in  obtaining  high-quality  solutions  to  (8),  using 
simulation  to  first  derive  a  convex  outer-bound  on  the  set  of  feasible  V  parameters  [17],  followed 
by  coordinatewise  optimization  over  V  and  (si, s25, 53<$)-  Nevertheless,  the  nonconvexity  is  not 
to  be  taken  lightly,  and  any  numerical  attempt  to  compute  pD{Vpoiy,S)  must  itself  treated  as  a 
lower  bound. 

Finally,  note  that,  if  l\  and  /2  have  positive  definite  quadratic  part,  then  the  feasibility  of 
(8)  implies  the  robust  stability  of  the  uncertain  linearized  dynamics  using  common  quadratic 
Lyapunov  function.  For  systems  with  cubic  vector  fields,  the  feasibility  of  (8)  is  also  necessary 
by  the  following  theorem  whose  proof  is  in  the  Appendix. 

Theorem  2.1:  Let  /0, . . . ,  fm  be  cubic  polynomials  in  x  satisfying  /o(0)  =  . . .  =  /m(0)  =  0, 
P  >~  0,  L\  >-  0,  L2  >-  0,  p(x)  =  xTPx,  l\(x )  =  xTL\x ,  and  l2(x)  =  xTL2x.  For  <5  e  A,  let 
A$  be  such  that  A$x  is  the  linear  (in  x)  part  of  fo(x)  +  YaL\^iMx)-  If  there  exists  Q  y  0 
satisfying  AjQ  +  QAs  -<  0  for  all  8  €  £a,  then  the  constraints  in  (8)  are  feasible.  o 

III.  POLYNOMIAL  PARAMETRIC  UNCERTAINTY 

We  extend  the  development  in  section  II  to  systems  with  polynomial  parametric  uncertainty. 
Specifically,  we  consider  the  system 


(9) 


771  TTLpu 

*(*) = /o(*«) + E  + E  »(*)  w*«). 

t=l  j-1 

where  /0,  /i,  -  -  -  7  /mj  fm+\  *  •  •  j  fm+mpu  :  lZn  — ►  lZn  are  vector  valued  polynomial  functions 
satisfying  /o(0)  =  . . .  =  /m+mpu(0)  =  0,  and  pi, €  R[<J]  are  scalar  valued  polynomial 
functions,  and  S  takes  values  in  a  bounded  poly  tope  A.  Note  that  pi, . . . ,  gmjni  are  bounded  since 
they  are  polynomials  with  the  bounded  domain  A.  We  begin  with  =  1  (for  simplicity)  and 
then  generalize  for  >  1. 

Replacing  pi  (6)  by  an  artificial  parameter  cp,  the  dynamics  in  (9)  can  be  written  as 

771 

*(0  =  fo(x(t))  +  '}T6ifi(x(t))  +  <t>fm+1(x(t)).  (10) 

»=i 

Our  approach  is  based  on  covering  the  graph  of  p,  {(£,  p(C))  €  7Zm+l  :  £  G  A}  ,  by  a  bounded 

polytope  T  C  7Zm+l.  Then,  the  dependence  of  the  vector  field  in  (10)  on  the  parameters  (5,<p) 

is  affine  and  ( 5 ,  <f> )  takes  values  in  the  bounded  polytope  I\  Therefore,  results  from  section  II 

are  applicable  for  the  system  in  (10)  by  replacing  A  by  T. 

A  polytope  T  covering  the  graph  of  p  can  be  obtained  by  bounding  p  from  above  and 

below  by  affine  functions  a%5  +  bu  and  aj5  +  bi  over  the  set  A,  namely  T (ai,au,bi,bu)  := 

{(£,  xp)  G  lZm+l  :  £  G  A,  af£  +  bi<xp<  a^£  +  6U}  .  The  volume  of  T  is  a  linear  function  of 

ai,au,bi ,  and  bu,  Volume(r(o;,  ou,  fy,  bu))  =  ( au  —  ai)T  /  £d£  4-  ( bu  -  bt)  /  The  polytope 

J  A  J  A 

with  smallest  volume  among  such  covering  polytopes  can  be  characterized  via 


min  Volume(r(a/,  au.  6/,  bu))  subject  to 

CL l 


(ID 


g(S)-(afS  +  bt)>  0  and  g(5)  -  (c%5  +  bu)  <  0,  V5  G  A. 

Using  the  generalized  S-procedure  [16],  an  upper  bound  for  this  minimal  volume  can  be  com¬ 
puted  by  a  linear  SOS  optimization  problem.  To  this  end,  let  affine  functions  hi,  i  =  1, ...  ,N, 
provide  an  inequality  description  for  A,  i.e.,  A  =  {£  G  7£m  :  /ii(C)  >0,  i  =  1, . . . ,  N}  . 
Proposition  3.1:  The  value  of  the  optimization  problem 

Volume(r(a/,  au,  bi,  bu ))  subject  to  aU{  G  £[<J],  an  G  E[5]  ,  i  =  1  ,...,N, 

-  g(8)  +  (al5  +  bu)  -  ZL  <7ui(6)hi(6)  G  E[<f],  (12a) 

9(6)  -  (af8  +  bt)  -  M6)hi(6)  G  E[5]  (12b) 


mm 

Ql,ClUyblybxi  ^ui^SuiyGueSli 


is  an  upper  bound  for  (1 1).  Here  <S’s  are  finite  dimensional  subsets  of  R[<J]. 


< 


Remarks  3.1:  Note  that  Volume(r(a/,au,6i,6u))  =  Volume(T(0,au,  0, 6U))-Volume(r(0,  a*,  0, 6*)), 
and  therefore  the  optimizing  values  of  the  variables  a*,au,  bi,  and  bu  in  Proposition  3.1  can 
equivalently  be  computed  by  two  smaller  optimization  problems.  < 

In  case  >  1,  affine  upper  and  lower  bounds  for  git . . . ,  gmpn  can  be  used  to  construct  a 
polytope  (with  2m+mPu  vertices)  covering  the  graph  of  (pi, . . . , gmpu ) ,  [1], 

Proposition  3.2:  For  j  =  1 - ,  m^,  let  a^S  +  by  and  a^8  +  buj  be  affine  functions  bound¬ 

ing  <7j  over  A  from  below  and  above,  respectively.  Then,  the  polytope  F  with  the  vertex 

set  £t  :=  (J  {(c,  ^1,  ■  ■  • , ^rnpu)  €  Km+mr™  :  +  bQj,  a  €  {/,  u},  j  =  1, . . . ,  m^} 

contains  the  graph  of  (pi , . . . ,  pmpu ) .  < 

This  gives  one  specific  procedure  to  cover  the  graph  of  a  vector-valued  multivariate  polynomial 
by  a  convex  polytope.  Further  research  that  advances  graph  covering  strategies  and  quantifies 
the  trade-off  between  the  number  of  vertices  and  the  volume  of  the  covering  polytope  would  be 
relevant  and  applicable  to  the  robust  ROA  problem.  Finally,  Proposition  3.2  can  be  used  with 
bounded  non-polynomial  p's  as  long  as  affine  upper  and  lower  bounds  are  provided. 

IV.  BRANCH-AND-BOUND  TYPE  REFINEMENT  IN  THE  PARAMETER  SPACE 

The  optimization  problem  in  (8),  when  applied  with  D  =  A,  provides  a  method  for  computing 
invariant  subsets  of  the  robust  ROA  characterized  by  a  single  Lyapunov  function.  Therefore, 
results  by  (8)  may  be  conservative:  the  certified  invariant  subset  may  be  too  small  relative  to  the 
robust  ROA.  On  the  other  hand,  a  less  conservative  estimate  of  the  robust  ROA  can  be  obtained 
by  solving  (8)  for  each  8  G  A  with  D  =  {5}.  For  a  subset  DC  A,  define 

Pb(Vpoiy,S)  :=  min  f3{s}(Vpoiy,S).  (13) 

q£D 

Then,  /3A(Vpoty,S)  <  P*^(Vpoiy,S).  However,  computing  /^(V^y,  S)  'requires  solving  an  opti¬ 
mization  problem  for  each  8  e  A,  and  consequently  is  impractical.  Next,  we  propose  an  informal 
“branch-and-bound”  type  procedure  for  computing  lower  and  upper  bounds  for  ^(V^,  S ),  i.e., 
localizing  the  value  of  /3*^(Vpoiy,S).  The  method  is  based  on  computing  a  different  Lyapunov 
function  for  each  cell  of  a  finite  partition,  V ,  of  A. 

Branch-and-bound  (B&B)  is  an  algorithmic  method  for  global  optimization  based  on  two 
steps:  first  the  search  region  is  partitioned  into  a  union  of  smaller  regions  ,  or  cells  (branching) 
and  then  upper  and  lower  bounds  for  the  objective  function  restricted  to  each  cell  are  computed 


(bounding)  [10].  These  steps  are  repeated,  refining  the  partition  each  repetition  (e.g.  subdividing 
the  cell  with  the  worst  lower  bound).  If  the  upper  and  lower  bounds  are  such  that  their  difference 
converges  to  zero  uniformly  as  the  size  of  the  cell  goes  to  zero,  then  the  B&B  algorithm  converges 
to  a  global  optimum.  Without  such  specific  guarantees,  steps  are  simply  repeated  until  the  gap 
between  the  upper  and  lower  bounds  gets  suitably  small  or  a  maximum  number  of  steps  is 
reached.  Additionally,  for  our  problem,  we  take  into  account  the  polytopic  covering  described 
in  section  HI,  and  recompute  this  covering  whenever  any  cell  is  subdivided. 

The  lower  and  upper  bounds  are  defined  over  any  polytope  D  E  A.  Certainly  Pd  (Vpo/y,«S) 
is  a  lower  bound  for  P*D(Vpoty,S)-  Upper  bounds  for  Pb(Vpoiy,S)  can  be  obtained  via  divergent 
trajectories  and  infeasibility  of  certain  necessary  conditions  for  the  constraints  in  (8).  Let  6  E  D 
and  pnc(S)  be  the  minimum  value  of  p  attained  on  all  non-convergent  trajectories  of  (5),  with 
Pnc(8)  :=  oo  if  there  is  no  non-convergent  trajectory.  Since  every  trajectory  entering  an  invariant 
subset  of  the  robust  ROA  has  to  converge  to  the  origin,  ilPypnc^  cannot  be  a  subset  of  the  robust 
ROA;  hence,  for  any  Vpoiy  and  5,  Pb(Vpoiy,S)  <  Pnc(5).  Unfortunately,  pnc( 8 ),  as  defined,  is 
impossible  to  compute.  But,  any  non-convergent  trajectory  yields  an  upper  bound  on  p^c{8),  and 
consequently  on  Pb(Vpoiy,  S).  In  order  to  establish  another  upper  bound,  let  P  >  0  and  8  €  D 
be  fixed.  If  there  exists  V  E  V  certifying  that  S2Pi(g  is  in  the  robust  ROA  through  the  constraints 
in  (7),  then  V  has  to  be  (i)  positive  for  all  nonzero  x  E  7 Zn,  (ii)  less  than  or  equal  to  1  (without 
loss  of  generality)  and  decreasing  along  every  trajectory  of  (5)  (for  this  fixed  8)  starting  in  (2Pi/3. 
Therefore,  if  no  V  E  V  satisfies  properties  (i)  and  (ii),  then  there  is  no  V  E  V  certifying  that 
ilp'P  is  in  the  robust  ROA  via  (7).  The  minimum  such  value,  denoted  plp(8),  is  an  upper  bound 
on  P*D{Vpoiy,S).  In  the  case  V  is  parameterized  as  V(x)  =  aTz(x)  with  z  a  vector  of  basis 
functions  and  a  a  vector  of  real  scalar  decision  variables,  constraints  on  V  along  trajectories  are 
affine  constraints  on  a;  consequently,  an  upper  bound  on  Plp(5)  can  be  determined  by  simulation 
and  linear  programming  (see  [17]  for  more  details).  As  in  all  B&B  algorithms,  the  minimum 
(over  the  subsets  that  make  up  the  partition  of  A)  of  these  upper  and  lower  bounds  are  upper 
and  lower  bounds  for  Pbi^poiyiS)- 

V.  IMPLEMENTATION  ISSUES 

The  optimization  problem  in  (8)  provides  a  recipe  to  compute  invariant  subsets  of  the  robust 
ROA.  However,  the  number  of  constraints  in  (8)  and  consequently  the  number  of  decision 


variables  increase  exponentially  with  m  +  because  (8d)  contains  a  SOS  constraint  for 
each  vertex  value  of  the  uncertainty  polytope.  The  increase  in  the  problem  size  may  render 
(8)  computationally  challenging  for  even  modest  values  of  m  +  m^.  This  difficulty  can  be 
partially  alleviated  by  accepting  suboptimal  solutions  for  (8)  in  a  sequential  manner  [16].  To  this 
end,  let  D  be  a  polytopic  subset  of  A  and  Dsampie  be  a  finite  sample  in  D: 

•  Solve  (8)  with  Dsampie ,  Vpoiy,  and  S  and  call  the  optimizing  Lyapunov  function  Vsampie. 

•  For  each  8  €  £d,  compute 

7 s  ■—  max  7  subject  to  s2s  €  Eta],  and  S3S  €  Eta], 

0<7.s2i€<S2,s3aeS3  ( 14^ 

—  [(7  —  Vsample)s2S  +  VVjamp|e(/o  +  8ifi))s3j  +  (2]  €  E[x], 
and  define  Yubopt  :=  min  {7,5  :  8  G  £d}  ■  At  this  point,  ^v,ampl  bopt  is  an  invariant  subset 
of  the  robust  ROA. 

•  Determine  the  largest  sublevel  set  £tp  ^«.6opt^Vpo(  s  of  p  contained  in  ^vJomp,  .i7^t,opt  by 
solving 

l3s£bopt (Vpoiy,  S)  :=  max  (3  subject  to  Si  €  £[x] 

sieslt0  (15) 

-to?  -  p)si  +  Vsample  -  Yubopt]  e  £[*]. 

While  this  sequential  procedure  sacrifices  optimality  (i.e.,  for  a  given  polytopic  subset  DC  A, 
,jsubopt^y^oiy^  ^  <  pD(Vpoiy,  <£)),  it  provides  practical  advantages:  For  a  fixed  Lyapunov  function 
candidate  Vsampie,  constraints  in  (8d)  (which  contain  one  SOS  constraint  for  each  vertex  value 
of  D )  decouple.  Therefore,  it  is  possible  to  determine  largest  value  of  7  such  that  ^v,ampU,y  C 
{x  €  lZn  :  ^Vsampie{x)f{x,8)  <  0}  for  every  8  €  £d  by  solving  (14)  independendy  for  each 
8  €  £d- 

Remarks  5.1:  Let  D  C  A  and  Dsampie  cDbea  singleton.  Then,  the  value  pDsarnpieC^poiy,  £>)— 
psvbopt^y^^  5)  js  always  non-negative  and  can  be  interpreted  as  a  measure  of  potential  improve¬ 
ment  in  the  lower  bound  for  ^(Vpoiy.S)  from  further  sub-dividing  D  in  the  DhB  refinement 
procedure.  Therefore,  it  may  be  used  as  a  stopping  criterion  in  an  informal  B&B  algorithm. 
However,  we  re-emphasize  that  pDt,ampU(ypoiy ,<S)  is  computed  solving  a  non-convex  optimization 
problem,  so  that  its  use  as  an  upper  bound  is  ad  hoc  and  referred  to  as  a  “quasi-upper”  bound 
(for  example  see  Figure  2).  <3 


Fig.  1.  Top  figures:  Bounds  for  /3f0,i]  vs*  number  of  B&B  iterations  with  d(V)  =  2  (left)  and  d(V)  =  4  (right).  Curves 
with  44o”  are  for  the  lower  bounds  obtained  by  directly  solving  (8)  with  D  taken  as  the  vertices  of  the  corresponding  cell  and 
curves  with  44o”  are  for  the  lower  bounds  obtained  by  applying  the  sequential  procedure  from  section  V  by  taking  Daamj>u  as 
the  center  of  the  corresponding  cell.  Bottom  figure:  Intersections  of  sublevel  sets  of  Vys  certified  to.be  in  the  robust  ROA  with 
d(V)  =  2  (inner  red,  solid  curve)  and  d(V)  =  4  (outer  red,  solid  curve).  Sublevel  sets  of  p  certified  to  be  in  the  robust  ROA 
d{V)  =  2  (inner  black,  dashed  curve)  and  d{V)  =  4  (outer  black,  dashed  curve),  estimate  of  the  robust  ROA  reported  in  [6] 
(blue,  dotted  curve).  Gray  dots  are  the  initial  conditions  of  trajectories  which  do  not  converge  to  the  origin  for  some  <5  6  (0, 1J. 


VI.  Examples 

For  the  following  examples,  we  implemented  the  sequential  procedure  from  section  V  using 
the  method  from  [17]  in  the  first  step  with  h(x)  =  l2(x)  =  10_6xrx  and  p(x)  =  xTx. 

A.  An  example  from  the  literature 

Consider  the  system,  [6],  governed  by 


X\ 

-6x2  +  x\  +  x\ 

5  + 

4x2  -  A 

r2 

X  = 

+ 

(T 

3x\  —  2x2 

-10xi  +  6x2  + 

12xi  —  4x2 

with  5  €  [0, 1].  We  applied  the  refinement  procedure  with  the  initial  partition  {[0, 1]}  for  d{V)  = 
2  and  d{V )  =  4.  Upper  and  lower  bounds  for  /3(*0  l,  (top  left  for  d{V)  =  2  and  top  right  for 
d(V)  =  4)  and  certified  invariant  subsets  of  the  robust  ROA  are  shown  in  Fig.  1.  In  both  cases, 


the  first  iteration  (a  parameter  independent  Lyapunov  function  for  A,  [16])  and  even  a  few  more 
do  not  yield  a  certified  region. 


B.  Controlled  short  period  aircraft  dynamics 

We  apply  the  robust  ROA  analysis  for  uncertain  controlled  short  period  aircraft  dynamics  (see 
the  Appendix  for  the  parameters  used  in  the  model) 


Coi(xp)  +  S^nixp)  +  S\q3  i(xp) 

^ bxp  +  ^11  +  612^1 

Xp  = 

902 (£p)  "b  ^1^12^?  ~b  S2y22 (Xp) 

+ 

621  +  b22S2 

Xl 

0 

where  xp  =  [x\  x2  x3]r ,  x\,  x2,  and  x3  denote  the  pitch  rate,  the  angle  of  attack,  and  the  pitch 
angle,  respectively,  coi  and  cn  are  cubic  polynomials,  <702,  <722,  and  y31  are  quadratic  polynomials, 
£\2  and  £b  are  vectors  in  723,  611,612,621,  and  622  G  72,  and  u,  the  elevator  deflection,  is  the 
control  input.  Variations  in  the  center  of  gravity  in  the  longitudinal  direction  are  modeled  by 
<5i  €  [0.99,2.05]  and  variations  in  the  mass  are  modeled  S2  G  [—0.1, 0.1]..  The  control  input  is 
determined  by  x4  =  — 0.864yi  H — 0.321y2  and  u  =  2x4,  where  x4  is  the  controller  state  and 
the  plant  output  y  =  [xi  x3]T .  Define  x  :=  [xp  x4] T .  We  applied  the  branch-and-bound  type 
refinement  procedure  with  d(V)  =  2  and  d(V )  =  4  using  the  sequential  implementation  on  a 
computer  cluster  with  9  processors:  after  the  first  BhB  iteration,  the  cell  with  the  smallest  lower 
bound  is  subdivided  into  3  subcells  and  cells  with  2-nd,  3-rd,  and  4-th  smallest  lower  bounds 
are  sub-divided  into  2  subcells.  Fig.  2  shows  the  lower  bounds  and  upper  bounds.  Smallest  value 
of  p  attained  on  divergent  trajectories,  /3nc,  is  8.60  and  obtained  for  (<5i,  S2)  =  (2.039,  —0.099) 
and  the  initial  condition  (0.17,2.65,  —0.10, 1.24). 


C.  Controlled  short  period  aircraft  dynamics  with  first-order  unmodeled  dynamics 

Consider  the  closed  loop  system  in  Figure  3  where  uncertain  first-order  dynamics  are  intro¬ 
duced  between  the  controller  output  ( v )  and  the  plant  input  (u)  from  section  VI-B 


u(s )  =  1.25  +  G(s;  S3,  <54)u(s)  = 


1.25  +  0.75<53 


s  +  5  a 


v{s). 


(16) 


Here,  S3  G  [—1,1]  and  S4  G  [10  2,102]  are  uncertain  parameters  and  G(s;  <53,  <54)  is  introduced 
to  examine  the  effect  of  unmodeled  dynamics  on  the  ROA.  Let  X5  =  —S4Xs  —  S4v  and  u  = 


Vv  "quasi-upper"  and  lower  bounds  for  d(V)  =  4 
A 


number  of  B&B  steps 

Fig.  2.  Lower  bounds  for  0^  with  d(V)  =  2  (solid  black  with  “x”)  and  d(V)  =  4  (solid  blue  curve  with  “o”)  and  0nc  (solid 
red  with  “o")  computed  at  the  centers  of  the  cells  generated  by  the  B&B  Algorithm  for  the  d(V)  ~  4  run.  Dashed  curves 
are  for  (computed  values  of)  0{sy  where  5  is  the  center  of  the  cell  with  the  smallest  lower  bound  at  the  corresponding  step  of 
the  BSzB  refinement  procedure  for  d(V)  =  2  (dashed  black  with  “x”)  and  d(V)  =  4  (dashed  blue  with  “o”). 


Fig.  3.  Closed-loop  system  with  the  uncertain  first-order  dynamics  between  the  controller  and  the  plant  ( Sp  =  (<5i ,  ^2))- 

1. 5*53X5  +  (1.25  +  0. 75<53)i>  be  a  realization  of  G  and  x  =  [xj  x\  X5] T  denote  the  state  of 
the  closed  loop  dynamics.  We  applied  the  BhB  refinement  with  d(V)  =  2  for  two  cases:  (i) 
(<5i,<52)  =  (1.52,0),  (ii)  <5i  and  <52  are  treated  as  uncertain  parameters  (as  in  section  VI-B)  where 
the  resultant  vector  field  is  affine  in  <5i,  *52,  <53,  <54,  <5i<53,  *52<53,  and  <5,  and  the  covering  polytopes 
are  in  7 Z7  with  128  vertices.  For  p(x)  =  xTx,  £^4  90  is  shown  to  be  in  the  robust  ROA  for  case 
(i)  whereas  £2P>2  so  is  certified  to  be  in  the  robust  ROA  for  case  (ii). 

VII.  CONCLUSIONS 

This  paper  considers  the  problem  of  finding  certified,  inner-estimates  of  the  region-of-attraction 
for  a  certain  class  of  uncertain  nonlinear  systems.  At  its  core,  the  solution  approach  combines 
Lyapunov  analysis,  S-procedure  relaxations,  and  SOS/SDP  optimization.  Four  factors  contribute 
to  the  problem  complexity:  number  of  state  variables;  degree  of  vector  field;  number  of  uncertain 
parameters;  dependence  of  vector  field  on  uncertain  parameters.  The  challenges  associated  with 
state  dimension  and  vector  field  degree  (often  large  optimization  problems)  appear  somewhat 


common  across  solution  techniques.  By  contrast,  the  issues  which  arise  from  uncertainty  are 
attacked  using  a  variety  of  diverse  techniques. 

We  address  the  difficulties  due  to  parameter  uncertainty  through  parallelization,  partitioning 
the  parameter  space,  solving  a  large  number  of  (uncoupled)  sub-problems.  While  the  Lyapunov 
function  for  each  sub-problem  is  independent  of  the  uncertain  parameter,  the  net  result  yields  a 
parameter-dependent  (piecewise-constant  in  the  parameter)  Lyapunov  function.  This  is  an  alter¬ 
native  to  more  direct  approaches  which  use  explicitly  parameter-dependent  Lyapunov  functions, 
e.g.  [13],  [5],  and  a  single  optimization  (with  additional  indeterminate  and  decision  variables, 
used  to  represent  the  uncertain  parameters  and  capture  their  constraints)  to  solve  the  problem. 

Of  course,  the  question  of  how  fine  the  parameter  space  partition  must  be  before  the  proposed 
method  yields  a  certified  robust  ROA  is  still  largely  open,  so  it  is  impossible  to  say  that  one 
approach  is  superior/inferior  to  another.  Similarly,  we  do  not  claim  that  the  proposed  strategy 
is  practical  for  all  instances  of  systems  modeled  by  (9).  Indeed,  large  numbers  of  uncertain 
parameters,  entering  the  dynamics  in  complex  ways  might  require  an  untenable  level  of  parameter 
space  partitioning  to  yield  a  positive  result.  Nevertheless,  we  have  illustrated  the  approach  on 
several  academic,  but  nontrivial,  examples,  including  a  5-state,  4-parameter  model  with  non-affine 
parameter  dependence.  Moreover,  for  cubic  (in  state)  vector  fields,  we  have  a  (weak)  positive 
result  which  follows  from  Theorem  2.1:  for  any  specific  partition  of  the  parameter  space,  if 
over  each  cell,  the  linearized  uncertain  dynamics  are  quadratically  stable,  then  the  certification 
conditions  (8)  are  guaranteed  to  be  feasible  (with  analytically  derived  choices  for  the  decision 
variables).  Among  other  things,  this  implies  that  the  uncertain  linearization  could  provide  insight 
into  the  level  of  parameter  space  division  needed  for  robust  region-of-attraction  certification. 
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IX.  Appendix 


Let  z(x )  be  a  vector  of  all  monomials  of  degree  2  with  no  repetition  and  nx  be  its  length. 
Lemma  9.1:  Let  Q  €  lZr  and  Q  =  QT  ^  0.  Then  there  exists  a  positive  definite  matrix 
H  €  1Zn*xnz  such  that  xTxxTQx  =  z(x)THz(x).  < 

Proof:  Let  Li  e  1Znxnz  be  such  that  L{z(x)  =  XiX ,  then  ( xTx )  xTQx  =  ^"=1(£ix)TQ(xix)  = 
S"=i  z(x)LfQLiz(x )  =  z(x)THz(x)  Note  that  L  —  [Lj . . .  L^JT  has  full  column  rank  since 
every  entry  of  z(x)  is  XjXk  for  some  1  <  j,  k  <  n;  consequently,  H  =  LT(In  <g>  Q)L  >■  0.  ■ 

Proof  of  Theorem  2.1:  Let  Q  +  0  satisfy  ASQ  +  QAS  +  —2 L2,  for  all  6  e  £A,  and  Q  h  Li 
(such  Q  can  be  obtained  by  properly  scaling  Q).  Let  e  =  Amjn(L2),  V(x)  xTQx ,  and  H  be  a 
positive  definite  Gram  matrix  for  ( xTx)V(x )  (which  exists  by  Lemma  9.1).  Let  M2s  €  lZnxn‘, 
and  the  symmetric  matrix  Mzs  €  lZnzXUz  be  such  that  xTM2sz(x),  and  z{x)T M^z{x)  are 
cubic  and  quartic  (in  x)  parts  of  W(/o(x)  +  ^7™. }  5t/»(x)),  respectively.  Define  Si(x)  = 
Amox((3)/Amm(^>),  =  Oi2SxT x  with  a2( 5  =  Amo*  (M3^  +  ±MfsM2S)  / Knin  ( H )  (where  for 

a  symmetric  matrix  A,  A+  denotes  the  projection  on  the  positive  semidefinite  cone),  s3(5(x)  =  1, 

7  =  min  {e/ (2q2s)  :  5  G  £A}  ,  and  f3  =  7/(2si).  Then,  V  -  li  and  -  [(/3  -  p)si  +  (V  -  7)] 
are  SOS  since  they  are  positive  semidefinite  quadratic  polynomials.  For  <5  €  £A,  bs(x)  = 
-  [(7  -  +  i£fsS36  +  h]  =  [xT  z{x)T]  Bs  [xT  z(x)T]T ,  where 


--ia25I-L2-(AjQ  +  QAe) 

-Mv/2 

f/  -M2S/ 2 

-Mll/2 

csl2sP1  —  Mzs 

_  -M?s/2  «2SH-Mzs_ 

(17) 

Note  that  ct2SH  =  H  h  Xmax  (M+  +  l(MfsM2S)  I  +  A ^  (Mm  +  I 

Mzs  +  i;MfsM2S.  Consequently,  Bs  is  positive  semidefinite  by  the  Schur  complement  formula 
applied  to  the  far  left  term  in  (17)  and  bs  £  £[x].  ■ 

Parameters  for  the  uncertain  controlled  short  period  aircraft  dynamics:  coi  (xp)  =  —0.24366x2+ 
0.082272xix2  +  0.30492x^  +  0.015426x2x3  -  3.1883X!  -  2.7258x2  -  0.59781x3;  £b  =  [0  - 

0.041136  0]T;6n  =  1.594150;  q02(xp)  =  -0.054444x^+0.10889x2x3-0.054444xl+0.91136x1- 
0.64516x2  -  0.016621x3;  621  =  0.0443215;  Cn(xp)  =  0.30765xi  +  0.099232x^  +  0.12404XJ  + 
0.90912x2 +  0.023258x3;  bu  =  -0.06202;  in  =  [0  0.00045754  0]T;  <?22(xp)  =  -0.054444x^  + 
0.10889x2x3  -  0.054444x3  -  0.6445x2  -  0.016621x3;  b22  =  0.044321. 
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Local  Stability  Analysis  for  Uncertain  Nonlinear  Systems 

Ufuk  Topcu  and  Andrew  Packard 

Abstract — We  propose  a  method  to  compute  provably  invariant  subsets 
of  the  region-of-attraction  for  the  asymptotically  stable  equilibrium 
points  of  uncertain  nonlinear  dynamical  systems.  We  consider  polynomial 
dynamics  with  perturbations  that  either  obey  local  polynomial  bounds 
or  are  described  by  uncertain  parameters  multiplying  polynomial  terms 
in  the  vector  field.  This  uncertainty  description  is  motivated  by  both 
incapabilities  in  modeling,  as  well  as  bilinearity  and  dimension  of  the 
sum-of-squares  programming  problems  whose  solutions  provide  invariant 
subsets  of  the  region-of-attraction.  We  demonstrate  the  method  on  three 
examples  from  the  literature  and  a  controlled  short  period  aircraft 
dynamics  example. 

1.  Introduction 

We  consider  the  problem  of  computing  invariant  subsets  of  the 
region-of-attraction  (ROA)  for  uncertain  systems  with  polynomial 
nominal  vector  fields  and  local  polynomial  uncertainty  description. 
Since  computing  the  exact  ROA,  even  for  systems  with  known 
dynamics,  is  hard,  researchers  have  focused  on  finding  Lyapunov 
functions  whose  sublevel  sets  provide  invariant  subsets  of  the  ROA 
[1],  [2],  [3],  [4],  [5].  Recent  advances  in  polynomial  optimization 
based  on  sum-of-squares  (SOS)  relaxations  [6],  [7]  are  utilized  to 
determine  invariant  subsets  of  the  ROA  for  systems  with  known 
polynomial  and/or  rational  dynamics  solving  optimization  problems 
with  matrix  inequality  constraints  [8],  [9],  [10],  [11],  [12],  [13].  Ref. 
[14]  provides  a  generalization  of  Zubov’s  method  to  uncertain  sys¬ 
tems  and  [15]  investigates  robustness  of  the  ROA  under  time- varying 
perturbations  and  proposes  an  iterative  algorithm  that  asymptotically 
gives  the  robust  ROA.  Parametric  uncertainties  are  considered  in  [16], 
[17],  [18].  The  focus  in  [16]  is  on  computing  the  largest  sublevel  set 
of  a  given  Lyapunov  function  that  can  be  certified  to  be  an  invariant 
subset  of  the  ROA.  [17],  [18]  propose  parameter-dependent  Lyapunov 
functions  which  lead  to  potentially  less  conservative  results  at  the 
expense  of  increased  computational  complexity. 

Similar  to  other  problems  in  local  analysis  of  dynamical  systems 
based  on  Lyapunov  arguments  and  SOS  relaxations  [9],  [11],  [12], 
[17],  [13],  [19],  our  formulation  leads  to  optimization  problems  with 
bilinear  matrix  inequality  (BM1)  constraints.  BMIs  are  nonconvex 
and  bilinear  SDPs  (those  with  BM1  constraints)  are  generally  harder 
than  linear  SDPs.  Consequently,  approaches  for  solving  SDPs  with 
BMIs  are  limited  to  local  search  schemes  [20],  [21],  [22],  [23]). 

The  uncertainty  description  detailed  in  section  III  contains  two 
types  of  uncertainty:  uncertain  components  in  the  vector  field  that 
obey  local  polynomial  bounds  and/or  uncertain  parameters  appearing 
affinely  and  multiplying  polynomial  terms.  Using  this  description, 
we  develop  an  SDP  with  BMIs  to  compute  robustly  invariant  subsets 
of  the  ROA.  The  number  of  BMIs  (and  consequently  the  number  of 
variables)  in  this  problem  increases  exponentially  with  the  sum  of  the 
number  of  components  of  the  vector  field  containing  uncertainty  with 
polynomial  bounds  and  the  number  of  uncertain  parameters.  One  way 
to  deal  with  this  difficulty  is  first  to  compute  a  Lyapunov  function  for 
a  particular  system  (imposing  extra  robustness  constraints)  and  then 
determine  the  largest  sublevel  set  in  which  the  computed  Lyapunov 
function  serves  as  a  local  stability  certificate  for  the  whole  family  of 
systems.  Once  a  Lyapunov  function  is  determined  for  the  system  in 
the  first  step,  second  step  involves  solving  smaller  decoupled  linear 
SDPs.  Therefore,  this  two  step  procedure  is  well  suited  for  parallel 
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computation  leading  to  relatively  efficient  numerical  implementation. 
Moreover,  recently  developed  methods  [13],  [24],  which  use  simula¬ 
tion  to  aid  in  the  nonconvex  search  for  Lyapunov  functions,  extend 
easily  to  the  robust  ROA  analysis  using  simulation  data  for  finitely 
many  systems  from  the  family  of  possible  systems  (e.g.  systems 
corresponding  to  the  vertices  of  the  uncertainty  polytope  when  the 
uncertainty  can  be  described  by  a  polytope).  In  the  examples  in  this 
paper,  we  implement  this  generalization  of  the  simulation  based  ROA 
analysis  method  from  [13],  [24]. 

The  rest  of  the  paper  is  organized  as  follows:  Section  11  reviews 
results  on  computing  the  ROA  for  systems  with  known  polynomial 
dynamics.  Section  111  is  devoted  to  the  discussion  of  the  motivation 
for  this  work  and  the  setup  for  the  uncertain  system  analysis.  In 
section  IV  provides  a  generalization  of  the  results  from  section 
11  to  the  case  of  dynamics  with  uncertainty.  The  methodology  is 
demonstrated  with  three  small  examples  from  the  literature  and  a 
five-state  example  in  section  V. 

Notation:  For  x  €  lZn ,  x  >z  0  means  that  Xk  >  0  for  k  =  1,  •  ■  • ,  n. 
For  Q  =  Qt  €  ttnxn,  Q  >z  0  (Q  >-  0)  means  that  xTQx  >  0 
(<  0)  for  all  x  €  7£n.  R[x]  represents  the  set  of  polynomials  in 
x  with  real  coefficients.  The  subset  £[x]  :=  {7r  €  R[x]  :  n  = 

fli  4-  tt2  H - f* 7Tm»  tti,  •  •  • ,  7rm  €  R[x]}  of  R[x]  is  the  set  of  SOS 

polynomials.  For  n  €  R[x],  <9(7r)  denotes  the  degree  of  tt.  For  subsets 
X\  and  X2  of  a  vector  space  A',  X\  -j-  X2  :=  {xi  +  xa  :  X\  € 
X\y  X2  €  X2}.  In  several  places,  a  relationship  between  an  algebraic 
condition  on  some  real  variables  and  state  properties  of  a  dynamical 
system  is  claimed,  and  same  symbol  for  a  particular  real  variable  in 
the  algebraic  statement  as  well  as  the  state  of  the  dynamical  system 
is  used.  This  could  be  a  source  of  confusion,  so  care  on  the  reader’s 
part  is  required.  < 

II.  Computation  of  invariant  subsets  of 
REGION-OF-ATTRACTION 

In  this  section,  we  give  a  characterization  of  invariant  subsets  of 
ROA  using  Lyapunov  functions  and  formulate  a  bilinear  optimization 
problem  for  computing  these  functions  when  they  are  restricted  to 
be  polynomial.  These  results  will  be  modified  to  compute  invariant 
subsets  of  the  ROA  for  systems  with  uncertainty  in  section  IV.  Now, 
consider  the  system  governed  by 

x(t)  =  /(x(t)),  (1) 

where  x(t)  €  7ln  is  the  state  vector  and  /  :  TV1  — ♦  7in  is  such 
that  /( 0)  =  0,  i.e.,  the  origin  is  an  equilibrium  point  of  (1)  and  /  is 
locally  Lipschitz  on  Tin.  Let  <p(t;xo)  denote  the  solution  to  (1)  with 
the  initial  condition  x(0)  =  xo-  If  the  origin  is  asymptotically  stable 
but  not  globally  attractive,  one  often  wants  to  know  which  trajectories 
converge  to  the  origin  as  time  approaches  00.  This  gives  rise  to  the 
following  definition  of  the  region-of-attraction: 

Definition  2.1:  The  region-of-attraction  Ro  of  the  origin  for  the 
system  (1)  is 

:=  (xo  e  7Zn  :  lim  (p(t\  xq)  =  0)  . 

I  t— *  00  ) 

< 

For  77  >  0  and  a  function  V  :  1Zn  — >  7 Z,  define  the  77-sublevel  set 
of  V  as 

fiv*  :=  {x  €  7ln  :  V(x)  <  77}. 

For  simplicity,  Hvu  is  denoted  by  Civ.  Lemma  2.2  provides  a 
characterization  of  invariant  subsets  of  the  ROA  in  terms  of  sublevel 
sets  of  appropriate  Lyapunov  functions. 

Lemma  2.2:  If  there  exists  a  continuously  differentiable  function 
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V  :  7 Zn  — ♦  7?.  such  that 

V(0)  =  0  and  V(x)  >  0  for  all  x  ^  0,  (2) 

Qv  is  bounded,  and  (3) 

Qv\{0}c{xeUn  :  VV(x)/(x)  <  0}  ,  (4) 

then  for  all  xo  G  Qv,  the  solution  of  (1)  exists,  satisfies  <p(t;x o)  G 
Hu  for  all  t  >  0,  and  limt— «>  <£>(£; xo)  =  0,  i.e.,  Qv  is  an  invariant 
subset  of  Ro.  < 

Lemma  2.2  is  proven  in  [11],  [12]  using  a  similar  result  from  [25]. 
If  the  dynamical  system  has  an  exponentially  stable  linearization,  one 
can  impose  a  stricter  condition  replacing  (4),  for  fi  >  0,  by 

Qv\  {0}  C  {x  €  7^n  :  VV(x)/(x)  <  -fiV{x)}  (5) 

With  nonzero  /z,  (5)  not  only  assures  that  trajectories  starting  in  Qv 
stay  in  Qv  and  converge  to  the  origin  but  also  imposes  a  bound 
on  the  rate  of  exponential  decay  of  V  certifying  the  convergence 
and  provides  an  implicit  threshold  for  the  level  of  a  disturbance  that 
could  drive  the  system  out  of  fiy.  Therefore,  one  may  consider  the 
stability  property  implied  by  (5)  with  nonzero  fi  to  be  more  desirable 
in  practice.  With  this  in  mind,  all  subsequent  derivations  contain  the 
fiV  term.  The  relaxed  condition  in  equation  (4)  can  be  recovered  by 
setting  /z  =  0. 

III.  Setup  and  Motivation 

We  now  introduce  the  uncertainty  description  used  in  the  rest  of  the 
paper  and  explain  its  usefulness  in  ROA  analysis  based  on  computing 
Lyapunov  functions  using  SOS  programming.  Consider  the  system 
governed  by 

x(£)  =  f(x(t))  =  fo(x(t))  +  <t>(x(t))  +  0(s(O),  (6) 

where  /o,0,  0  :  7 Zn  — >  7Zn  are  locally  Lipschitz.  Assume  that  /0  is 
known,  0  €  £><*> ,  and  0  G  where 

:=  {0  :  (x)  <  0(x)  <  0u(x)  Vx  €  £}, 

V ^  :=  {0  :  0(x)  =  ^(x)q  Vx  €  G,  oti  ^  a  d:  <^u}- 

Here,  Q  is  a  given  subset  of  7 Zn  containing  the  origin,  0z  and  0U  are 
n  dimensional  vectors  of  known  polynomials  satisfying  0z(x)  X  0  < 
<t>u{x)  for  all  x  €  £,  a,az,au  €  7£N,  and  ^  is  a  matrix  of  known 
polynomials.  Let  0i,  0z,»,  0ti,t,  ai,t,  and  au,i  denote  i-th  entry 
of  0,  0i,  0u,  a,  c*z,  and  au,  respectively.  Define  V  :=  £>*-(- We 
assume  that  /o(0)  =  0,  0(0)  =  0  for  all  0  G  £><*>  (i.eM  0z(O)  =  0, 
and  0^(0)  =  0),  and  0(0)  =  0  for  all  0  G  P*,  (i.e.,  tf(0)  =  0)» 
which  assures  that  all  systems  in  (6)  have  a  common  equilibrium 
point  at  the  origin.1  In  order  to  be  able  to  use  SOS  programming,  we 
restrict  our  attention  to  the  case  where  /o,  0z,  4>u>  and  ^  have  only 
polynomial  entries  and  Q  is  defined  as  Q  :=  {x  G  7£n  :  g(x)  >z 
0,  gi  €  R[x],  i  =  1, . . .  ,m}.  Note  that  entries  of  0  do  not  have  to 
be  polynomial  but  have  to  satisfy  local  polynomial  bounds. 

Motivation  for  this  kind  of  system  description  stems  from  the 
following  sources: 

(i)  Perturbations  as  in  (6)  may  be  due  to  modeling  errors,  aging, 
disturbances  and  uncertainties  due  to  environment  which  may 
be  present  in  any  realistic  problem.  Prior  knowledge  about 
the  system  may  provide  local  bounds  on  the  entries  of  0 
and/or  bounds  for  the  parametric  uncertainties  a.  Moreover, 

*The  assumption  thai  all  possible  systems  in  (6)  have  a  common  equilibrium 
poini  can  be  alleviaied  by  generalizing  the  analysis  based  on  contraction 
metrics  and  SOS  programming  studied  in  [26]  to  address  local  stability 
(raiher  than  global  stability  as  in  [26]).  However,  lhis  method  leads  to  higher 
computaiional  cosl.  Therefore,  we  do  not  pursue  this  direction  here. 


TABLE  I 

NsdP  (LEFT  COLUMNS)  AND  Ndecision  (RIGHT  COLUMNS)  FOR 
DIFFERENT  VALUES  OF  n  AND  2d. 


2d 

n 

4 

6 

8 

10 

2 

6 

6 

10 

27 

15 

75  - 

21  165 

5 

21 

105 

56 

1134 

126 

6714 

'  252  2e4 

9 

55 

825 

220 

le4 

715 

2e5 

★  ★ 

14 

120 

4200 

680 

★ 

★ 

★ 

★  ★ 

16 

153 

6936 

i  ★ 

★ 

★ 

★ 

★  ★ 

uncertainties  that  do  not  change  system  order  can  always  be 
represented  as  in  (6)  (see  p.339  in  [27]). 

(ii)  Analysis  of  dynamical  systems  using  SOS  programming  is 
often  limited  to  systems  with  polynomial  or  rational  vector 
field.  In  [28],  a  procedure  for  re-casting  non-rational  vector 
fields  into  rational  ones  at  the  expense  of  increasing  the  state 
dimension  is  proposed.  Another  way  to  deal  with  a  non¬ 
polynomial  vector  field  is  to  locally  approximate  the  vector 
field  with  a  polynomial,  and  bound  the  error.  For  practical 
purposes  only  finite  number  of  terms  can  be  used.  Finite-term 
approximations  are  relatively  accurate  in  a  restricted  region 
containing  the  origin.  However,  they  are  not  exact.  On  the  other 
hand,  it  may  be  possible  to  represent  terms,  for  which  the  error 
between  the  exact  vector  field  and  its  finite-term  approximation 
obey  local  polynomial  bounds,  using  0  in  (6). 

(iii)  SOS  programming  can  be  used  to  analyze  systems  with 
polynomial  vector  fields.  The  number  of  decision  variables 
^decision  and  the  size  Ns  dp  of  the  matrix  in  the  SDP  for 
checking  existence  of  a  SOS  decomposition  for  a  degree  2d 
polynomial  in  n  variables  grows  polynomially  with  n  if  d  is 
fixed  and  vice  versa  [6].  However,  /Vsdp  and  /Vdecbion  get 
practically  intractable  for  the  state-of-the-art  SDP  solvers  even 
for  moderate  values  of  n  for  fixed  d  (see  Table  2,  where  solid 
lines  in  the  table  represent  a  fuzzy  boundary  between  tractable 
and  intractable  SDPs).  Moreover,  using  higher  degree  Lyapunov 
functions  and/or  higher  degree  multipliers  (used  in  the  sufficient 
conditions  for  certain  set  containment  constraints  in  section  IV) 
as  well  as  higher  degree  vector  fields  increases  the  problem 
size,  and,  in  fact,  the  growth  of  the  problem  size  with  the 
simultaneous  increase  in  n  and  d  is  exponential.  Therefore, 
in  order  to  be  able  to  use  SOS  programming,  one  may  have 
to  simplify  the  dynamics  by  truncating  higher  degree  terms  in 
the  vector  field.  In  this  case,  0z  and  0U  provide  local  bounds 
on  the  truncated  terms.  This  is  discussed  further  at  the  end 
of  section  (IV).  It  is  also  worth  mentioning  that  bilinearity, 
a  common  feature  of  the  optimization  problems  for  local 
analysis  using  Lyapunov  arguments  (see  section  IV),  introduces 
extra  complexity  [29]  and  therefore  a  further  necessity  for 
simplifying  the  system  dynamics. 

In  summary,  representation  in  (6)  and  definitions  of  V $  and  T)^ 
are  motivated  by  uncertainties  introduced  due  to  incapabilities  in 
modeling  and/or  analysis. 

IV.  Computation  of  robustly  invariant  sets 

In  this  section,  we  will  develop  tools  for  computing  invariant 
subsets  of  the  robust  ROA.  The  robust  ROA  is  the  intersection  of  the 
ROAs  for  all  possible  systems  governed  by  (6)  and  formally  defined, 
assuming  that  the  origin  is  an  asymptotically  equilibrium  point  of  (6) 
for  all  $  G  D,  as 
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Definition  4.1.  The  robust  ROA  of  the  origin  for  systems 
governed  by  (6)  is 

Bq  :=  Pi  {xo  G  Rn  :  lim  y{t\ xo, 5)  =  0}, 

1  1  £ — *  oo 

<5ez> 

where  <p(*;xo,<5)  denotes  the  solution  of  (6)  for  S  £  V  with  the 
initial  condition  x(0)  =  xo.  < 

The  robust  ROA  is  an  open  and  connected  subset  of  Rn  containing 
the  origin  and  is  invariant  under  the  flow  of  all  possible  systems 
described  by  (6)  [1 5].  We  focus  on  computing  invariant  subsets  of  the 
robust  ROA  characterized  by  sublevel  sets  of  appropriate  Lyapunov 
functions.  Since  the  uncertainty  description  for  <p  and  holds  only 
for  x  £  £,  we  will  also  require  the  computed  invariant  set  to  be  a 
subset  of  Q.  To  this  end,  we  modify  Lemma  2.2  such  that  condition 
(5)  holds  for  (6)  for  all  S  £  V  (i.e.,  for  all  <f>  £  T>t  and  rp  €  T>^). 

Proposition  4.2:  If  there  exists  a  continuously  differentiable  func¬ 
tion  V  :  7 Zn  — ►  7 Z  and  p  >  0  such  that,  for  all  6  £  V,  conditions 

(2M3). 

Qv  C  Q,  and  (7) 

Qv\ {0}  C  {x  €  7^n  :  VV(x)(f0(x)  +  S(x))  <  -pV(x)}  (8) 

hold,  then  for  all  x0  £  Hv  and  for  all  5  £  Vy  the  solution 
of  (6)  exists,  satisfies  <p(*;xo,<5)  £  Hv  for  all  t  >  0,  and 
lim£_oo  <p(t\ xo,  6)  =  0,  i.e.,  f lv  is  an  invariant  subset  of  Rq.  < 
Proof:  Proposition  4.2  follows  from  Lemma  2.2.  Indeed,  for 
any  given  system  x  =  /o(x)  4*  <5(x),  (8)  assures  that  (5)  is  satisfied. 
Then,  for  any  fixed  S  £  V  and  for  all  xo  £  Ov,  <p(£;xo,<5)  exists 
and  satisfies  <p(t;  xo,  <5)  £  Hv  for  all*  >  0,  limf—oo  xo,  (5)  =  0, 
and  Qv  is  an  invariant  subset  of  {xo  £  Rn  :  <p(t\xoy6)  — *  0}. 
Therefore,  Qv  is  an  invariant  subset  of  Rq.  ■ 

Remark  4.3:  In  fact,  Qv  is  invariant  for  both  time- invariant  and 
time-varying  perturbations.  The  conclusion  of  Proposition  4.2  holds 
for  time-varying  <5(=  4>  4-  a)  as  long  as  4>i(x)  <  <t>(xyt)  ■<  fa(x) 
and  ai  ■<  a(t)  ■<  au  for  all  x  £  Q  and  t  >  0.  Recall  that,  in  the 
uncertain  linear  system  literature,  e.g.  [30],  the  notion  of  quadratic 
stability  is  similar,  where  a  single  quadratic  Lyapunov  function  proves 
the  stability  of  an  entire  family  of  uncertain  linear  systems.  < 
Note  that  V  has  infinitely  many  elements;  therefore,  there  are 
infinitely  many  constraints  in  (8).  Now,  define 

St  ;=  {<p  :  fa  £  R[x]  and  fa  is  equal  to  fa%i  or  fa,i}y 
£ t  :=  {ip  :  ip{x)  =  ^(x)q,  ch  is  equal  to  or  au,i}, 


(3),  (7),  (9),  and  V0  C  Qv.  This  can  be  written  as  an  optimization 
problem. 

0*  (V)  max  0  subject  to 

V  €V,/3>0  J 

1/(0)  —  0  and  V(x)  >  0  for  all  x  /  0, 

Qv  =  {x  £  RP  :  V{x)  <  1 }  is  bounded 
V0  =  {x£Rn  :  p(x)  <0}C  QVy 
Qv  eg, 

nv\{  o>c 

n  ix  6  :  V'/(/oW  +  £(*))  <  “M* 

Here,  V  denotes  the  set  of  candidate  Lyapunov  functions  over  which 
the  maximum  is  defined  (e.g.  V  may  be  equal  to  all  continuously 
differentiable  functions). 

In  order  to  make  the  problem  in  (10)  amenable  to  numerical 
optimization  (specifically  SOS  programming),  we  restrict  V  to  be 
a  polynomial  in  x  of  fixed  degree.  We  use  the  well-known  sufficient 
condition  for  polynomial  positivity  [6]:  for  any  tt  £  R[x],  if  n  £ 
E[x],  then  7r  is  positive  semidefinite.  Using  simple  generalizations  of 
the  5-procedure  (Lemmas  8.1  and  8.2  in  the  appendix),  we  obtain 
sufficient  conditions  for  set  containment  constraints.  Specifically,  let 
li  and  1%  be  a  positive  definite  polynomials  (typically  exTx  with 
some  (small)  real  number  e).  Then,  since  li  is  radially  unbounded, 
the  constraint 

V  —  l\  £  E[x]  (11) 

and  1/(0)  =0  are  sufficient  conditions  for  the  constraints  in  (10b). 
By  Lemma  8.1,  if  sj  £  E[x]  and  s^k  £  E[x]  for  k  =  1 .  ,m,  then 

-[(0-p)sx  +  (V-l)]€£[x]  (12) 

gk  -(1  -  V)s4fc  £  £[x],  k =  1, . . . ,  m,  (13) 

imply  the  first  and  second  constraints  in  (10c),  respectively.  By 
Lemma  8.2,  if  S2(,s$(  £  E[x]  for  £  £  £,  then 

-  [(1  -  V)s2t  4  (VV(f0  4  0  4  pV)sH  4-  h]  £  E[x]  (14) 

is  a  sufficient  condition  for  the  feasibility  of  the  constraint  in  (lOd). 
Using  these  sufficient  conditions,  a  lower  bound  on  0*(V)  can  be 
defined  as  an  optimization: 

Proposition  4.5:  Let  0B  be  defined  as 

0B{Vpoiy,S)  :=  max  0  subject  to  (11)  —  (14),  (15) 


l, 


(10a) 

(10b) 

(10c) 


'(*)}.  00d) 


and  S  :=  St  4-  S^,.  5,  a  finite  subset  of  D,  can  be  used  to  transform 
condition  (8)  to  a  finite  set  of  constraints  that  are  more  suitable  for 
numerical  verification: 

Proposition  4.4:  If 

Qv\  {0}  C  {x  £  7^n  :  W(x)(/0(x)  +  t(x))  <  -pV(x)}  (9) 

holds  for  all  ^  £  5,  then  (8)  holds  for  all  5  £  V.  < 

Proof:  Let  x  £  Qv  be  nonzero  and  5  £  T>.  Then,  x  £  Q  by  (7); 
therefore,  there  exist  l\ ,  •  •  • ,  tn,  k\ , . . . ,  fcjv  (depending  on  x)  with 
0  <  ti  <  1  and  0  <  ki  <  1  such  that  <5(x)  =  Lfa(x)  4*  (I  — 
L)fa(x)  4-  V{x){Kai  4 -  (/  -  if)au),  where  L  and  K  are  diagonal 
with  La  =  £i  and  Ku  =  ki.  Hence,  there  exist  nonnegative  numbers 
i/(  (determined  from  £'s  and  fc’s)  for  (  £  8  with  Ylfze  =  ^ 
such  that  (5(x)  =  Consequently,  by  VU(x)(/0(x)  4 

S(x))  =  W(f)(/0(x)+E€€f  =  E<ec^V(*)(/ b(*)+ 

C(5))  <  (8)  follows.  ■ 

In  order  to  enlarge  the  computed  invariant  subset  of  the 
robust  ROA,  we  define  a  variable  sized  region  V0  := 
{x  £  RP  :  p(x)  <  /3},  where  p  £  R[x]  is  a  fixed,  positive  definite, 
convex  polynomial,  and  maximize  0  while  imposing  constraints  (2)- 


V(0)  =  0,  V  £  Vpoly ♦  3l  £  5l,  3  €  520  S34  £  53^, 
S4fc  6  54fc,  and  0  >  0.  Here,  Vpoiv  C  V  and  5’s  are  prescribed 
finite-dimensional  subspaces  of  R[x]  and  E[x],  respectively  Then, 
0B(Vpoly,S)<0*(Vpoly).  < 

The  optimization  problem  in  (15)  provides  a  recipe  to  compute 
subsets  of  Rn  that  are  invariant  under  the  flow  of  all  possible  systems 
described  by  (6).  The  number  of  constraints  in  (15)  (consequendy  the 
number  of  decision  variables  since  each  new  constraint  includes  new 
variables)  increases  exponentially  with  N  and  n  —  n  where  n  is 
defined  as  the  number  of  entries  of  the  vectors  <t>i  and  fa.  satisfying 
fa(x)  —  fa(x)  =  0  for  all  x  £  Q.  Namely,  there  are  2(n~n)+N 
SOS  conditions  in  (15)  due  to  the  constraint  in  (14).  Revisiting  the 
discussion  in  item  (iii)  at  the  end  of  section  Ill,  we  note  that  covering 
the  high  degree  vector  field  with  low  degree  uncertainty  reduces 
the  dimension  of  the  SOS  constraints  but  increases  (exponentially, 
depending  on  n  -  n)  the  number  of  constraints.  Consequently,  the 
utility  of  this  approach  will  depend  on  n  —  n  and  is  problem 
dependent.  Example  (3)  in  section  V-A  illustrates  this  technique. 

This  difficulty  can  be  partially  alleviated  by  accepting  suboptimal 
solutions  for  (15)  in  two  steps:  First  compute  a  Lyapunov  function  for 
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a  finite  sample  of  systems  corresponding  to  the  finite  set  V sample  C 
V  (for  example,  T>8ampU  can  be  taken  as  the  singleton  corresponding 
to  the  “center”  of  V)  solving  the  problem 


max  0  subject  to 
V-h€  E[ar], 

'  -K0-P)*i  +  (V—  1)]€E[*]  <16> 

9k  -  (1  -  V)s4*  €  E(x],  fc= 

-  ((1  -  V)s2*  +  W(/o  +  6)s36  +  la]  €  £[x], 

for  S  €  'P sample  ^  where  3i  €  *Si,S26  €  <^2,335  €  <$3,  €  *^4 

are  SOS,  V  €  Vpoiv>  V’(O)  =  0,  and  let  V  the  Lyapunov  function 
computed  by  solving  (16).  In  the  second  step,  compute  the  largest 
sublevel  set  V0»uboPt  such  that  V  certifies  V0»»b*Pt  to  be  in  the  ROA 
for  every  vertex  system  by  solving  several  smaller  decoupled  affine 
SDPs.  For  £  €  £,  define 

7^  :=  max  7  subject  to  32,33  €  £[x], 

7.«3€52,«3€53 

-  [(7  -  V)s2  +  VV(/o  +  0«3  + 12]  €  E[x], 

and  Yubopt  min{7^  :  £  €  £}.  Then,  a  lower  bound  for  08Ub°vt 
can  be  computed  through 


psubopt  max  p  subject  to  si  €  E[x] 

€S  1 


-  [(/3-p)Sl  +  (V-7*“t>0''t)]  €E[x]. 


(18) 


While  the  two-step  procedure  sacrifices  optimality,  it  has  prac¬ 
tical  computational  advantages.  The  constraints  in  (14)  decouple 
in  the  problem  (17).  In  fact,  for  each  £  €  £a,  the  problem  in 
(17)  contains  only  a  single  constraint  from  (14).  Therefore,  this 
decoupling  enables  suboptimal  local  stability  analysis  for  systems 
with  uncertainty  without  solving  optimization  problems  larger  than 
those  one  would  have  to  solve  in  local  stability  analysis  for  systems 
without  uncertainty.  Furthermore,  problems  in  (17)  can  be  solved 
independently  for  different  f  G  £a  and  therefore  computations  can 
be  trivially  parallelized.  Advantages  of  this  decoupling  may  be  better 
appreciated  by  noting  that  one  of  the  main  difficulties  in  solving 
large-scale  SDPs  is  the  memory  requirements  of  the  interior-point 
type  algorithms  [31].  Consequently,  it  is  possible  to  perform  some 
ROA  analysis  on  systems  with  relatively  reasonable  number  of  states 
and/or  uncertain  parameters  using  the  proposed  suboptimal  solution 
technique. 


Finally,  the  following  upper  bound  on  the  value  of  py  for  which 
(14)  can  be  feasible,  will  be  useful  in  section  V. 


Proposition  4.6.  Let  L2  >-  0  and  h(x)  :=  xT  L2X  Then, 

p.  := 


max  p.  subject  to 

M>0,P=Pr>;0 

AjP  +  PA{  +  9.P1  -L2,  for  all  £  €  £, 


(19) 


where  := 
such  that  (14)  can  be 


±Q 


,  is  an  upper  bound  for  the  values  of  p 
feasible.  < 


Proof:  With  I2  as  defined  and  S2$,53<  €  E[x],  if  £>e(x)  := 
-[(1-  V)87z  +  (W(f0  +  0  +  +  l*\  €  E[x],  then  s2* 

and  cannot  contain  constant  and  linear  monomials  and  the 
quadratic  part  of  b$  has  to  be  SOS  and  equivalently  positive  semidef- 
inite.  Therefore,  the  result  follows  from  the  fact  that,  for  fixed  p  >  0 
and  positive  definite  quadratic  I2,  the  existence  of  P  y  0  satisfying 
AjP  4-  PA{  -F  pP  -Z/2  is  necessary  for  the  existence  of  V,  32$, 
and  33$  feasible  for  (14).  ■ 

Note  that  the  problem  in  (19)  can  be  solved  as  a  sequence  of  affine 
SDPs  by  a  line  search  on  p. 


Fig.  1.  Invariant  subsets  of  ROA  reported  in  [16]  (solid)  and  those  compuied 
solving  the  problem  in  (15)  with  d(V)  =2  (dash)  and  d(V)  =  4  (dot)  along 
with  initial  conditions  (stars)  for  some  divergeni  trajectories  of  the  system 
corresponding  to  a  =  1. 


TABLE  II 

Optimal  values  of  0  in  the  problem  ( 1 5)  with  different  values 
OF  p  AND  d(V)  =  2  AND  4. 


2 

4 

0 

0.623 

0.771 

0.01 

0.603 

0.763 

0.05 

0.494 

0.742 

0.1 

0.404 

0.720 

0.15 

0.277 

0.698 

0.2 

0.137 

0.676 

v.  Examples 

In  the  following  examples,  p(x)  =  xTx  (except  for  example  (2) 
in  section  V-A),  /i(x)  =  10“  6xTx,  and  h(x)  =  10“6xTx.  All 
certifying  Lyapunov  functions  and  multipliers  are  available  at  [32]. 
All  computations  use  the  generalization  of  the  simulation  based  ROA 
analysis  method  from  [13],  [24].  Representative  computation  times 
on  2.0  GHz  desktop  PC  are  listed  with  each  example. 

A.  Examples  from  the  literature 

(1)  Consider  the  following  system  from  [16]:  xi  =  X2  and 
X2  =  —X2  4*  q(— xi  -F  xf),  where  a  €  [1,3]  is  a  parametric 
uncertainty.  We  solved  problem  (15)  with  d(V)  =  2  and  d{V)  =  4 
for  p  =  0,0.01,0.05,0.1,0.15,  and  0.2.  Note  that  p  (as  defined 
in  Proposition  4.6)  is  0.244.  Typical  computation  times  are  5  and  8 
seconds  for  d{V)  =  2  and  4,  respectively. 

Figure  1  shows  the  invariant  subset  of  the  robust  ROA  reported 
in  [16]  (solid)  and  those  computed  here  with  d(V)  =  2  (dash)  and 
d(F)  =  4  (dot)  for  p  =  0  along  with  two  points  (stars)  that  are  initial 
conditions  for  divergent  trajectories  of  the  system  corresponding  to 
a  =  1.  Table  II  shows  the  optimal  values  of  0  in  the  problem  (15) 
with  d{V)  =  2  and  4  for  different  values  of  p. 

(2)  Consider  the  system  (from  [17])  of  xi  =  — X2  -F  0.2ax2  and 
X2  =  x\  -F  (x?  —  l)x2  where  a  €  [—1,1].  For  easy  comparison  with 
the  results  in  [17],  let  p(x)  =  0.378xi  —  0.274xiX2  +  0.278x2  and 
p  =  0.  In  [17],  it  was  shown  that  Po.545  (with  a  single  parameter 
independent  quartic  V)y  V^.ni  (with  pointwise  maximum  of  two 
parameter  independent  quartic  F’s),  Vq.goo  (with  a  single  parameter 
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Fig.  3.  Invariant  subsets  of  ROA  with  d{V)  =  2  (solid)  and  8(V)  =  4 
(dash)  along  with  initial  conditions  for  divergent  trajectories  ( “*  ”  for  4>{x)  = 
(0.76i|,0  19(i2+I2))and  for <t>(x)  =  (-0.76x|, -0.19(i?+i|))). 


Fig.  2.  Invariant  subsets  of  ROA  with  d(V )  =  4  (inner  solid)  and  d(  V)  =  6 
(dash)  along  with  the  unstable  limit  cycle  (outer  solid  curves)  of  the  system 
corresponding  to  a  =  — 1.0,  —0.8,  . .  ,0.8, 1.0. 


TABLE  III 

Optimal  values  of  (3  in  the  problem  (15)  with  different  values 

OF  p  AND  d(V)  =  4  AND  6. 


4 

6 

0 

0.773 

0.826 

0.01 

0.767 

0.820 

0.05 

0.741 

0.803 

0.1 

0.708 

0.787 

0.2 

0.640 

0.750 

0.5 

0.517 

0.651 

0.75 

0.406 

0.573 

dependent  quartic  (in  state)  V ),  Po  soe  (with  pointwise  maximum  of 
two  parameter  dependent  quartic  (in  state)  1^’s)  are  contained  in  the 
robust  ROA.  On  the  other  hand,  the  solution  of  problem  (15)  with 
d(V)  =  4  and  d(V)  =  6  certifies  that  P0.773  and  Po.826  are  subsets 
of  the  robust  ROA,  respectively.  Figure  2  shows  invariant  subsets  of 
the  robust  ROA  computed  using  d(V)  =  4  (inner  solid)  and  d(V)  = 
6  (dash)  along  with  the  unstable  limit  cycle  (outer  solid  curves)  of  the 
system  corresponding  to  a  =  —1.0,  —0.8, . . .  ,0.8, 1.0.  In  order  to 
demonstrate  the  effect  of  the  parameter  p  on  the  size  of  the  invariant 
subsets  of  the  robust  ROA  verifiable  solving  the  optimization  problem 
in  (15),  the  analysis  is  repeated  with  p  =  0.01, 0.05, 0.1.,  0.2, 0.5, 
and  0.75.  Note  that  p  (as  defined  in  Proposition  4.6)  is  0.769.  Table 
III  shows  the  optimal  values  of  (3  in  the  problem  (15)  with  d(V)  =  4 
and  6  for  different  values  of  p.  Typical  computation  times  are  19  and 
24  seconds  for  d{V)  —  4  and  6,  respectively. 

(3)  Consider  the  system  governed  by 


— 2xi  4-  £2  4*  x?  -I-  1.58x2 
— xi  —  £2  4-  0.13x2  4*  0.66x?£2 


+  0(x),  (20) 


where  <p  satisfies  the  bounds  -0.76x1  <  <p i(x)  <  0.76x1  and 
-0.19(x?  +x!)  <  <t> 2(x)  <  0.19(x?  +xl)  in  the  set  Q  =  {x  € 
V?  :  g(x)  =  xTx  <  2.1}.  Fig.  3  shows  invariant  subsets  of 
the  robust  ROA  computed  with  d(V)  =  2  (solid)  and  d(V)  =  4 
(dash)  along  with  two  points  that  are  initial  conditions  for  divergent 
trajectories  (“*”  for  <f>(x)  =  (0.76x1, 0.19(x?  +x!))  and  “  x  ”  for 
<j>(x)  =  (-0.76x1,  -0,19(x?  +  x!))).  Typical  computation  times  are 


13  and  35  seconds  for  d(V)  =  2  and  4,  respectively. 


B.  Controlled  short  period  aircraft  dynamics 
Consider  the  plant  dynamics 


-3  -1.35  -0.56  " 

1.35-  0.0422  ' 

2  = 

0.91  -0.64  -0.02 

2  4- 

0.4 

1  0  0 

0 

(1  4*  ai)(0.082iz2  4-  0.44z|  +  O.OI2223  4-  0.222|) 

+ 

(l  4-  a2)(-0.0522  4-  O.II2223  -  O.O523) 

0 

y  =  [z\  23] T.  where  z\y  22  and,  23  are  the  pitch  rate,  the  angle  of 
attack,  and  the  pitch  angle,  respectively  The  input  u  is  the  elevator 
deflection  and  determined  by 


'  -0.60 

0.09  ' 

V  + 

'  -0.06 

-0.02  ' 

0 

0 

-0.75 

-0.28 

u  =  7/1  4-2.2772,  where  77  is  the  controller  state.  Here,  a\  and  Q2  are 
two  uncertain  parameters  introducing  10%  uncertainty  for  the  entries 
of  the  plant  dynamics  that  are  nonlinear  in  t\  i.e.,  Qi  €  [—0.1, 0.1] 
and  a 2  €  [—0.1, 0.1].  Entries  in  the  vector  fields  above  are  shown 
up  to  three  significant  digits.  The  exact  vector  field  used  for  this 
example  is  available  at  [32].  The  solution  of  (15)  with  ^(V^)  =  2 
and  p  =  0  verifies  that  V7.2  C  Ro  whereas  it  can  be  certified 
that  V8.6  is  a  subset  of  the  ROA  for  the  nominal  system  (i.e.,  for 
Qitl  =  Qf, 2  =  a u.i  =  c*u,2  =  0).  With  d(V)  =  4  the  problem  in 
(15)  has  more  than  4000  decision  variables.  Therefore,  we  computed 
a  suboptimal  solution  in  two  steps  for  p  —  0:  We  first  computed 
a  Lyapunov  function  for  the  nominal  system  (35  minutes,  which 
certifies  that  V\ 5.2  is  in  the  ROA  for  the  nominal  system)  and  then 
verified  (3  minutes)  that  V9  6  is  an  invariant  subset  of  the  ROA  for 
the  uncertain  system.  To  assess  the  suboptimality  of  the  results,  we 
performed  extensive  simulations  for  the  uncertain  system  setting  Qi 
and  c*2  to  their  limit  values  and  found  a  diverging  trajectory  with 
the  initial  condition  satisfying  p(2(0),  77(0))  »  14.  The  gap  between 
the  value  of  (3  «  14  for  which  V0  cannot  be  a  subset  of  the  robust 
ROA  and  the  value  of  (3  =  9  6  for  which  C  Ro  is  verified 
may  be  due  to  the  finite  dimensional  parametrization  for  V ,  the 
issues  mentioned  in  Remark  4.3,  the  fact  that  we  only  use  sufficient 
conditions  and/or  suboptimality  of  the  two  step  procedure  used  for 
this  example;  nevertheless,  it  demonstrates  a  necessity  of  further  study 
to  make  local  system  analysis  based  on  Lyapunov  functions  and  SOS 
relaxations  more  efficient. 


6 


VI.  Conclusions 

We  proposed  a  method  to  compute  provably  invariant  subsets  of  the 
region -of-attracti on  for  the  asymptotically  stable  equilibrium  points 
of  uncertain  nonlinear  dynamical  systems.  We  considered  polynomial 
dynamics  with  perturbations  that  either  obey  local  polynomial  bounds 
or  are  described  by  uncertain  parameters  multiplying  polynomial 
terms  in  the  vector  field.  This  uncertainty  description  is  motivated  by 
both  incapabilities  in  modeling,  as  well  as  bilinearity  and  dimension 
of  the  sum-of- squares  programming  problems  whose  solutions  pro¬ 
vide  invariant  subsets  of  the  region-of-attraction.  We  demonstrated 
the  method  on  three  examples  from  the  literature  and  a  controlled 
short  period  aircraft  dynamics  example. 
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VIII.  Appendix 

Following  two  lemmas  are  simple  generalizations  of  the  S- 
procedure.  The  proof  of  the  first  one  is  trivial.  We  provide  a  proof 
for  the  second  one. 

Lemma  8.1:  Given  9o,9\,"  -  ,9m  €  R[x],  if  there  exist 

6  E[x]  such  that  go  -  Y^iLis*9<  €  £[x],  then 
{x  €  7Zn  :  pi  (x),. ..  ,gm{x)  >  0}  C  {x  e  7ln  :  g0(x)  >  0}  .  < 

Lemma  8.2:  Let  g  €  R[x]  be  positive  definite,  h  €  R[x], 
7  >  0,  51,52  6  E[x],  l  6  R[x]  be  positive  definite  and  satisfy 
/( 0)  =  0.  Suppose  that  —  [(7  —  g)si  +  hs2  +  I]  G  £[x]  holds.  Then, 
&gn\{ 0}  C  {x  €  7Zn  :  h(x)  <  0  and  52  (x)  >  0}  < 

Proof:  Let  x  €  be  nonzero  Then, 

0  >  -l(x)  -  (7  -  g(x))s 2(x)  >  h(x)s2(x), 
and,  consequently,  52 (x)  >  0  (since  52 (x)  >  0)  and  h(x)  <0.  ■ 
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Abstract 

The  problem  of  computing  bounds  on  the  region-of-attraction  for  systems  with  polynomial  vector  fields  is  considered.  Invariant 
subsets  of  the  region-of-attraction  are  characterized  as  sublevel  sets  of  Lyapunov  functions.  Finite  dimensional  polynomial 
parameter izations  for  Lyapunov  functions  are  used.  A  methodology  utilizing  information  from  simulations  to  generate  Lyapunov 
function  candidates  satisfying  necessary  conditions  for  bilinear  constraints  is  proposed.  The  suitability  of  Lyapunov  function 
candidates  are  assessed  solving  linear  sum-of-squares  optimization  problems.  Qualified  candidates  are  used  to  compute  invariant 
subsets  of  the  region-of-attraction  and  to  initialize  various  bilinear  search  strategies  for  further  optimization.  We  illustrate  the 
method  on  small  examples  from  the  literature  and  several  control  oriented  systems. 
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1  Introduction 

The  region-of-attraction  (ROA)  of  a  locally  asymptot¬ 
ically  stable  equilibrium  point  is  an  invariant  set  such 
that  all  trajectories  emanating  from  points  in  this  set 
converge  to  the  equilibrium  point.  Computing  the  exact 
ROA  for  nonlinear  dynamics  is  very  hard  if  not  impossi¬ 
ble.  Therefore,  researchers  have  focused  on  determining 
invariant  subsets  of  the  ROA.  Among  all  other  meth¬ 
ods  those  based  on  Lyapunov  functions  are  dominant  in 
the  literature  (Davison  and  Kurak,  1971;  Genesio  et  ai, 
1985;  Vannelli  and  Vidyasagar,  1985;  Chiang  and  Thorp, 
1989;  Chesi  et  al .,  2005;  Papachristodoulou,  2005;  Tan 
and  Packard,  2006;  Hachicho  and  Tibken,  2002;  Tibken 
and  Fan,  2006;  Tibken,  2000).  These  methods  compute 
a  Lyapunov  function  as  a  local  stability  certificate  and 
sublevel  sets  of  this  Lyapunov  function,  in  which  the 
function  decreases  along  the  flow,  provide  invariant  sub¬ 
sets  of  the  ROA. 

Using  sum-of-squares  (SOS)  relaxations  for  polynomial 
nonnegativity  (Parrilo,  2003),  it  is  possible  to  seaxch  for 
polynomial  Lyapunov  functions  for  systems  with  poly¬ 
nomial  and/or  rational  dynamics  using  semidefinite  pro¬ 
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gramming  (Papachristodoulou,  2005;  Tan  and  Packard, 
2006;  Hachicho  and  Tibken,  2002).  Reliable  and  efficient 
solvers  for  linear  semidefinite  programs  (SDPs)  are  avail¬ 
able  (Sturm,  1999).  However,  the  SOS  relaxation  for  the 
problem  of  computing  invariant  subsets  of  the  ROA  leads 
to  bilinear  matrix  inequality  (BMI)  constraints.  BMIs 
are  nonconvex  and  bilinear  SDPs,  those  with  BMI  con¬ 
straints,  are  known  to  be  NP-hard  in  general  (Toker  and 
Ozbay,  1995).  Consequently,  the  state-of-the-art  of  the 
solvers  for  bilinear  SDPs  is  far  behind  that  for  the  linear 
ones.  Recently  PENBMI,  a  solver  for  bilinear  SDPs,  was 
introduced  (Kocvara  and  Stingl,  2005)  and  subsequently 
used  for  computing  invariant  subsets  of  the  ROA  (Tern 
and  Packard,  2006;  Tibken  and  Fan,  2006).  It  is  a  local 
optimizer  and  its  behavior  (speed  of  convergence,  qual¬ 
ity  of  the  local  optimal  point,  etc.)  depends  on  the  point 
from  which  the  optimization  starts. 

By  contrast,  simulating  a  nonlinear  system  of  moderate 
size,  except  those  governed  by  stiff  differential  equations, 
is  computationally  efficient.  Therefore,  extensive  simu¬ 
lation  is  a  tool  used  in  real  applications.  Although  the 
information  from  simulations  is  inconclusive,  i.e.,  cannot 
be  used  to  find  provably  invariant  subsets  of  the  ROA,  it 
provides  insight  into  the  system  behavior.  For  example, 
if,  using  Lyapunov  arguments,  a  function  certifies  that  a 
set  V  is  in  the  ROA,  then  that  function  must  be  positive 
and  decreasing  on  any  solution  trajectory  initiating  in  V. 
Using  a  finite  number  of  points  on  finitely  many  conver¬ 
gent  trajectories  and  a  linear  parametrization  of  the  Lya- 
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punov  function  U,  those  constraints  become  affine,  and 
the  feasible  polytope  (in  ^-coefficient  space)  is  a  convex 
outer  bound  on  the  set  of  coefficients  of  valid  Lyapunov 
functions.  It  is  intuitive  that  drawing  samples  from  this 
set  to  seed  the  bilinear  SDP  solvers  may  improve  the 
performance  of  the  solvers.  In  fact,  if  there  are  a  large 
number  of  simulation  trajectories,  samples  from  the  set 
often  are  suitable  Lyapunov  functions  (without  further 
optimization)  themselves.  Effectively,  we  are  relaxing  the 
bilinear  problem  (using  a  very  specific  system  theoretic 
interpretation  of  the  problem)  to  a  linear  problem,  and 
the  true  feasible  set  is  a  subset  of  the  linear  problem’s 
feasible  set.  Information  from  simulations  is  also  used  in 
(Prokhorov  and  Feldkamp,  1999)  and  (Serpen,  2005)  for 
computing  approximate  Lyapunov  functions. 

Notation:  For  x  £  7 Zn,  x  X  0  means  that  Xk  >  0  for 
k  =  1,---  , n.  For  Q  =  QT  £  7lnxn,  Q  X  0  (Q  y  0) 
means  that  xTQx  >  0  (<  0)  for  all  x  £  7ln.  R[x]  repre¬ 
sents  the  set  of  polynomials  in  x  with  real  coefficients. 

The  subset  E[x]  :=  {n  £  R[x]  :  tt  =  +  7r|  + - h 

TTm,  » •  •  *  ,  7rm  €  R[x] }  of  R[x]  is  the  set  of  SOS  polyno¬ 
mials.  For  7r  G  R[x],  d(n)  denotes  the  degree  of  n.  Cl  de¬ 
notes  the  space  of  continuously  differentiable  functions. 
We  use  the  term  “semidefinite  programming”  to  mean 
optimization  problems  with  affine  objective  function  and 
general  (not  necessarily  affine)  matrix  (semi) definiteness 
constraints.  < 

2  Characterization  of  invariant  subsets  of  the 
ROA  and  bilinear  SOS  problem 

Consider  the  autonomous  nonlinear  dynamical  system 

x(t)  =  /(*(<))»  (1) 

where  x(Z)  £  7Zn  is  the  state  vector  and  /  :  7Zn  — >  1Zn 
is  such  that  /( 0)  =  0,  i.e.,  the  origin  is  an  equilibrium 
point  of  (1),  and  /  is  locally  Lipschitz.  Let  <£(£,  t)  denote 
the  solution  to  (1)  at  time  t  with  the  initial  condition 
</>(£,  0)  =  £.  If  the  origin  is  asymptotically  stable  but  not 
globally  attractive,  one  often  wants  to  know  which  tra¬ 
jectories  converge  to  the  origin  as  time  approaches  oo. 
The  region-of- attraction  R$  of  the  origin  for  the  system 
(1)  is  Rq  :=  {£  G  lZn  :  lim^oo  <£(£,  t)  =  0} .  A  modifica¬ 
tion  of  a  similar  result  in  (Vidyasagar,  1993)  provides 
a  characterization  of  invariant  subsets  of  the  ROA.  For 
rj  >  0  and  a  function  V  :  7Zn  — ►  1Z ,  define  the  77-sublevel 
set  of  V  as  fiv/y  *=  €  7Zn  :  V’(x)  <  77}. 

Lemma  1  Let  7  G  72.  be  positive .  If  there  exists  a  C 1 
function  V  :  1Zn  — >  71  such  Z/iaZ 


Qyi7  w  bounded ,  and  (2) 

1^(6)  =  0  and  V(x)  >  0  /or  all  x  £  lZn  (3) 

nV|7\  {0}  c{xG^n  :  W(x)/(x)  <  0}  ,  (4) 

i/ien  /or  a//  £  G  solution  of  (1)  exists ,  satisfies 

<£(£,  t)  £  fly, 7  /or  aZZ  Z  >  0,  and  lim(_oo  <£(£,  £)  =  0,  i.e., 
fiy,7  is  an  invariant  subset  of  Rq.  < 


In  order  to  enlarge  the  computed  invariant  sub¬ 
set  of  the  ROA,  we  define  a  variable  sized  region 
Vp  :=  {x  G  lZn  :  p(x)  <  /?},  where  p  £  R[x]  is  a  fixed 
positive  definite  convex  polynomial,  and  maximize  0 
while  imposing  the  constraint  Vp  C  Qy n  along  with  the 
constraints  (2)- (4).  This  can  be  written  as 

0*  (V)  :=  ^  max  ^  0  subject  to  (2)  —  (4),  Vp  C 

(5) 

Here  V  denotes  the  set  of  candidate  Lyapunov  func¬ 
tions  over  which  the  maximum  is  defined,  for  example  all 
Cl  functions.  Lemma  1  and  the  associated  optimization 
problem  in  (5)  provide  a  characterization  of  the  invari¬ 
ant  subsets  of  the  ROA  in  terms  of  the  sublevel  sets  of 
Lyapunov  functions. 

The  problem  in  (5)  is  an  infinite  dimensional  problem. 
In  order  to  make  it  amenable  to  numerical  optimization 
(specifically  SOS  optimization),  we  restrict  V  to  be  all 
polynomials  of  some  fixed  degree.  We  use  the  well-known 
sufficient  condition:  for  any  n  £  R[x],  if  n  £  E[x],  then 
7 r  is  positive  semidefinite  (Parrilo,  2003).  Using  simple 
generalizations  of  the  S-procedure  (Lemmas  2  and  3), 
we  obtain  sufficient  conditions  for  set  containment  con¬ 
straints.  Specifically,  let  Zi  and  l2  be  a  positive  definite 
polynomials  (typically  exTx  for  some  small  real  number 
e).  Then,  since  l\  is  radially  unbounded,  the  constraint 

V  —  l\  £  E[x]  (6) 

and  U(0)  =  0  are  sufficient  conditions  for  (2)  and  (3). 
By  Lemma  2,  if  S\  £  E[x],  then 

-[09-  P)*i  +  (V-  7)]  €E[x]  (7) 

implies  the  set  containment  Vp  C  and  by  Lemma 
3,  if  s2,  S3  G  E[x],  then 

-  [(7  -  V)s2  +  W/s3  +  h)  €  E[x]  (8) 

is  a  sufficient  condition  for  (4).  Using  these  sufficient 
conditions,  a  lower  bound  on  0*(V)  can  be  defined  as 

0*r(V,S)  :=  max  0  subject  to  (6)  —  (8), 

"  7  V€V,0,st€«st  (9) 

V(0)  =  0,  Sj  G  E[x],  /?  >  0. 

Here,  the  sets  V  and  Si  are  prescribed  finite-dimensional 
subspaces  of  polynomials.  Although  0*B  depends  on  these 
subspaces,  it  will  not  always  be  explicitly  notated.  Note 
that  since  the  conditions  (6)- (8)  are  only  sufficient  con¬ 
ditions,  /?b(V,«S)  <  /?*(V)  <  0 *(Cl).  The  optimization 
problem  in  (9)  is  bilinear  because  of  the  product  terms 
0S\  in  (7)  and  V  s2  and  W fs 3  in  (8).  However,  the  prob¬ 
lem  has  more  struct  me  than  a  general  BMI  problem.  If 
V  is  fixed,  the  problem  becomes  affine  in  5  =  {s i,  s2 ,  S3} 
and  vice  versa.  In  section  3,  we  will  construct  a  convex 
outer  bound  on  the  set  of  feasible  V  and  sample  from 
this  outer  bound  set  to  obtain  candidate  V’s,  and  then 
solve  (9)  for  S,  holding  V  fixed. 


3  Relaxation  of  the  bilinear  SOS  problem  using 
simulation  data 

The  usefulness  of  simulation  in  understanding  the  ROA 
for  a  given  system  is  undeniable.  Faced  with  the  task  of 
performing  a  stability  analysis  (eg.,  “for  a  given  p,  is  Vp 
contained  in  the  ROA?”),  a  pragmatic,  fruitful  and  wise 
approach  begins  with  a  linearized  analysis  and  at  least  a 
modest  amount  of  simulation  runs.  Certainly,  just  one  di¬ 
vergent  trajectory  starting  in  Vp  certifies  that  Vp  £  Ro. 
Conversely,  a  large  collection  of  only  convergent  trajec¬ 
tories  hints  to  the  likelihood  that  indeed  Vp  C  Ro-  Sup¬ 
pose  this  latter  condition  is  true,  let  C  be  the  set  of  Nconv 
trajectories  c  converging  to  the  origin  with  initial  condi¬ 
tions  in  Vp.  In  the  course  of  simulation  runs,  divergent 
trajectories  d  whose  initial  conditions  are  not  in  Vp  may 
also  get  discovered,  so  let  the  set  of  d’s  be  denoted  by 
D  and  N<tiV  be  the  number  of  elements  of  D.  Although 
C  and  D  depend  on  /3  and  the  manner  in  which  Vp  is 
sampled,  this  is  not  explicitly  notated. 

With  (3  and  7  fixed,  the  set  of  Lyapunov  functions  which 
certify  that  Vp  C  Ro,  using  conditions  (6)- (8),  is  sim¬ 
ply  {V  G  R[x)  :  (6)  -  (8)  hold  for  some  s*  G  E[x]}  .  Of 
course,  this  set  could  be  empty,  but  it  must  be  contained 
in  the  convex  set  {V  G  R[x]  :  (10)  holds},  where 

W(c(*))/(c(t))  <  0, 

h  (C (*))  <  nc(0),  and  V(c(0))  <  7,  (10) 

7  +  <5  <  V(d(*)), 

for  all  c  G  C,  d  G  D,  and  t  >  0,  where  6  is  a  fixed  (small) 
positive  constant.  Informally,  these  conditions  simply  say 
that  any  V  which  verifies  that  Vp  C  Ro  using  the  condi¬ 
tions  (6)- (8)  must,  on  the  trajectories  starting  in  Vp ,  be 
decreasing  and  take  on  values  between  0  and  7.  More¬ 
over,  V  must  be  greater  than  7  on  divergent  trajecto¬ 
ries.  In  fact,  with  the  exception  of  the  strengthened  lower 
bound  on  V  (beyond  mere  positivity),  the  conditions  in 
(10)  are  even  necessary  conditions  for  any  V  G  C* l  which 
verify  Vp  C  Ro  using  conditions  (2)-(4). 

3.1  Affine  relaxation  using  simulation  data 

Let  V  be  linearly  parameterized  as  V  :=  {V  G 
R[x]  :  V{x)  =  p(x)Ta},  where  a  G  Rnb  and 
<p  is  n^-dimensional  vector  of  polynomials  in  x . 
Given  <p(x),  constraints  in  (10)  can  be  viewed  as 
constraints  on  a  G  7Znb  yielding  the  convex  set 
{a  G  7lnb  :  (10)  holds  for  V  =  tp(x)Ta}.  For  each 

c  G  C,d  G  D,  let  Tc  and  7d  be  finite  subsets  of  the 
interval  [0,oo)  including  the  origin.  A  polytopic  outer 
bound  for  this  set  described  by  finitely  many  constraints 
is  3^aim  {ck  G  1Znb  :  (11)  holds},  where 

[Vv?(c(tc))/(c(tc))]Tq!  <  0, 

Mc(tc))  <  <p(c(Tc))Ta,  and  </j(c(0))Ta:  <  7,  (11) 

<iC’(d(rd))Ta  >7  +  5 


for  all  c  €  C,  rc  e  Tc,  d  e  D,  and  rd  €  Td.  Note 
that  ip(c(0))Ta  <  7  in  (11)  provides  necessary  condi- 
tions  for  Vp  C  QyyJ  since  c(0)  G  Vp  for  all  c  G  C. 

In  practice,  we  replace  the  strict  inequality  in  (11)  by 
[Vcp(c(rc))/(c(rc))]T  a  <  —  h{c(Tc)),  where  I3  is  a  fixed, 
positive  definite  polynomial  imposing  a  bound  on  the 
rate  of  decay  of  V  along  the  trajectories. 

The  constraint  that  W f  be  negative  on  a  sublevel  set 
of  V  implies  that  Wf  is  negative  on  a  neighborhood 
of  the  origin.  While  a  large  number  of  sample  points 
from  the  trajectories  will  approximately  enforce  this, 
in  some  cases  (eg.  exponentially  stable  linearization)  it 
is  easy  to  analytically  express  as  a  constraint  on  the 
low  order  terms  of  the  polynomial  Lyapunov  function. 
For  instance,  assume  V  has  a  positive-definite  quadratic 
part,  and  that  separate  eigenvalue  analysis  has  estab¬ 
lished  that  the  linearization  of  (1)  at  the  origin,  i.e., 
x  —  V/(0)x,  is  asymptotically  stable.  Define  C(P) 

(V/(0))TP  +  P(V/(0)),  where  PT  =  P  y  0  is  such 
that  x"Px  is  the  quadratic  part  of  V.  Then,  if  (8)  holds, 
it  must  be  that 

C(P)  0.  (12) 

Let  yun  :=  {a  G  Vnb  :  P  =  PT  X  0  and  (12)  holds}. 

It  is  well-known  that  yun  is  convex  (Boyd  and  Vanden- 
berghe,  2004).  Again,  in  practice,  (12)  is  replaced  by  the 
condition  C(P)  •<  — el ,  for  some  small  real  number  e. 
Furthermore,  define  D^sos  :=  {ot  G  7lnb  :  (6)  holds}. 
By  (Parrilo,  2003),  ysos  is  convex.  Since  34m  and 
3^sos  are  convex,  y  :=  ysim  fl  fl  ysos  is  a  convex 
set  in  Vnb.  Equations  (11)  and  (12)  constitute  a  set  of 
necessary  conditions  for  (6)- (8);  thus,  we  have  y  D  B  := 

{a  G  7Znb  :  3s2,S3  G  E[x]  such  that  (6)  —  (8)  hold}. 
Since  (8)  is  not  jointly  convex  in  V  and  the  multipliers,  B 
may  not  be  a  convex  set  and  even  may  not  be  connected. 

A  point  in  y  can  be  computed  solving  an  affine  (feasi¬ 
bility)  SDP  with  the  constraints  (6),  (11)  and  (12).  An 
arbitrary  point  in  y  may  or  may  not  be  in  B.  However,  if 
we  generate  a  collection  A  of  Ny  points 

distributed  approximately  uniformly  in  y ,  it  may  be  that 
some  of  the  points  are  in  B.  To  this  end,  we  use  the  so- 
called  “Hit-and-Run”  (H&cR)  random  point  generation 
algorithm  as  described  in  (Tempo  et  al.,  2005).  When 
applied  to  generate  a  sample  of  3^,  each  step  of  H&cR 
algorithm  requires  solving  four  small  affine  SDPs. 

3.2  Algorithms 

Since  a  feasible  value  of  0  is  not  known  a  priori,  an  iter¬ 
ative  strategy  to  simulate  and  collect  convergent  and  di¬ 
vergent  trajectories  is  necessary.  This  process  when  cou¬ 
pled  with  the  H&R  algorithm  constitutes  the  Lyapunov 
function  candidate  generation. 

Simulation  and  Lyapunov  function  generation  (SimLFG) 

algorithm:  Given  positive  definite  convex  p  G  R[x],  a 
vector  of  polynomials  ip(x)  and  constants  A^conv, 

AV,  ft  shrink  £  (0, 1),  and  empty  sets  C  and  D,  set  7=1, 

A^rnore  =  A ^conv,  ^div  ~  0* 

i.  Integrate  (1)  from  Nmore  initial  conditions  in  the 


set  {x  G  1Zn  :  p(x)  =  &sim}- 

ii.  If  there  is  no  diverging  trajectory,  add  the  trajec¬ 
tories  to  C  and  go  to  (iii).  Otherwise,  add  the  di¬ 
vergent  trajectories  to  D  and  the  convergent  tra¬ 
jectories  to  C,  let  Nd  denote  the  number  of  diverg¬ 
ing  trajectories  found  in  the  last  run  of  (i)  and  set 
Ndiv  to  Ndiv  +  Nd •  Set  0sim  to  the  minimum  of 
fishrinkPsiM  and  the  minimum  value  of  p  along  the 
diverging  trajectories.  Set  Nmare  to  Nmore  -  Nd , 
and  go  to  (i). 

iii.  At  this  point  C  has  Nc(mv  elements.  For  each  i  — 
1,  •  •  • ,  Nconv,  let  t i  satisfy  c <(r)  G  VpSIM  for  all 
r  >  Ti .  Eliminate  times  in  %  that  are  less  than  r* . 

iv.  Find  a  feasible  point  for  (6),  (11),  and  (12).  If 
(6),  (11),  and  (12)  are  infeasible,  set  0sim  = 
Pshrink&siM,  and  go  to  (iii).  Otherwise,  go  to  (v). 

v.  Generate  Ny  Lyapunov  function  candidates  using 

H&R  algorithm,  and  return  0sim  and  Lyapunov 
function  candidates.  < 

The  suitability  of  a  Lyapunov  function  candidate  is  as¬ 
sessed  by  solving  two  optimization  problems.  Both  prob¬ 
lems  require  bisection  and  each  bisection  step  involves  a 
linear  SOS  problem.  Alternative  linear  formulations  ap¬ 
pear  in  the  appendix.  These  do  not  require  bisection,  but 
generally  involve  higher  degree  polynomial  expressions. 

Problem  1:  Given  V  G  R[x]  (from  SimLFG  algorithm) 
and  positive  definite  I2  G  R[x],  define 

7 r  :=  max  7  subject  to  $2,53  €  Sfxl,  7  >  0, 

7,52.53  (13) 

-  ((7  -  *0*2  +  VV7s3  +  h]  €  £[x]. 

If  Problem  1  is  feasible,  then  7^  >  0  and  define 
Problem  2:  Given  V  G  R[x],  p  G  R[x],  and  7^,  solve 

0*L  :=  max  0  subject  to  s\  G  £[x],  0  >  0, 

0,51  (14) 

Although  72  and  0*L  depend  on  the  allowable  degree  of 
s  1,  52,  and  53,  this  is  not  explicitly  notated. 

Assuming  Problem  1  is  feasible,  it  is  true  that  Vp*  \{0}  C 
\{0}  C  {x  G  lln  :  VVr(x)/(x)  <  0},  so  V  certifies 
that  Vp*L  C  Rq.  Solutions  to  Problems  1  and  .2  provide  a 
feasible  point  for  the  problem  in  (9).  This  feasible  point 
can  be  further  improved  by  solving  the  problem  in  (9) 
using  PENBMI  and/or  iterative  coordinate- wise  linear 
optimization  schemes,  one  of  which  is  given  next. 

Coordinate-wise  optimization  (CWOpt)  algorithm : 

Given  V  G  R[x],  positive  definite  l\,h  G  R[x],  a  con¬ 
stant  eiter,  and  maximum  number  of  iterations  Nxter, 
set  h  =  0 

i.  Solve  Problems  1  and  2. 

ii.  Given  s  1,  52,  53,  and  72  from  step  (i),  set  7  in  (7)- 

(.8)  to  72,  solve  (9)  for  V  and  0,  and  set  0*L  =  0B- 

iii.  If  k  =  Niter  or  the  increase  in  0*L  between  successive 


applications  of  (ii)  is  less  than  e<£cr,  return  V ,  72, 
and  0*L.  Otherwise,  set  k  to  k  +  1  and  go  to  (i).  < 

The  algorithms  ( SimLFG ,  Problems  1  and  5,  and 
CWOpt)  yield  lower  bounds  on  0*(C* l),  as  they  produce 
a  Lyapunov  function  which  certifies  that  a  particular 
value  of  0  satisfies  Vp  G  Rq.  Upper  bounds  (i.e.,  values 
of  0  that  are  not  certifiable)  may  also  be  obtained.  More 
specifically,  diverging  trajectories  found  in  the  course  of 
simulation  runs  provide  upper  bounds  on  0*(Cl)  while 
inconsistency  of  the  constraints  (6),  (11),  and  (12)  pro¬ 
vide  upper  bounds  on  0*B.  A  diverging  trajectory  with 
the  initial  condition  xo  satisfying  p(xo)  =  0  proves  that 
Vp  cannot  be  a  subset  of  the  ROA,  i.e.,  0*  (C1 *)  <  0.  Fur¬ 
thermore,  restricting  Lyapunov  function  candidates  to 
V<p  :=  {< p(x)Ta  :  a  G  1Znbj  has  additional  implications. 
Infeasibility  of  any  of  the  constraints  (6),  (11),  and  (12) 
for  some  value  of  0  (recall  (11)  implicitly  depends  on 
0)  verifies  0*B  (V^,5)  <  0 *  (V^)  <  0 ,  regardless  of  the 
subspaces  constituting  S.  Moreover,  the  gap  between 
the  value  of  0  proven  unachievable  and  what  we  actu¬ 
ally  certify,  namely  a  lower  bound  to  0*B  (V^,,  <S),  can  be 
used  as  a  measure  of  suboptimality  introduced  due  to 
the  finiteness  of  the  degree  of  the  multipliers  and  the 
fact  that  the  bilinear  search  and  the  coordinate- wise  lin¬ 
ear  search  are  only  local  optimization  schemes.  Finally, 
H&R,  SimLFG  and  CWOpt  algorithms  become  more 
efficient  using  parallel  computing. 

4  Examples 

Certifying  Lyapunov  functions,  multipliers  and  missing 
parameters  for  all  examples  in  this  paper  are  available 
at  http: //j agger .me .berkeley . edu/~pack/certif y. 
In  the  examples,  i,*(x)  =  10~6xrx  for  i  —  1, 2, 3. 

4.1  Van  der  Pol  dynamics 

The  Van  der  Pol  dynamics  xi  =  — X2,  X2  =  Xi  4-  (x^  — 
l)x2  have  a  stable  equilibrium  point  at  the  origin  and  an 
unstable  limit  cycle.  The  limit  cycle  is  the  boundary  of 
the  ROA.  We  applied  SimLFG  algorithm  with  p(x)  = 
xTx  and  the  parameters  Ncorxv  =  200,  0sim  =30  (ini¬ 
tial  value),  0shrink  =  0.9,  and  Ny  —  50  for  d(V)  =  2,4, 
and  6.  We  found  Ndiv  =  21  diverging  trajectories  dur¬ 
ing  the  simulation  runs  and  feasible  solutions  for  (6), 
(11),  and  (12)  in  step  (iv)  with  0sim  =  1.44,  1.97, 
and  2.19  for  d(V)  =  2,4,  and  6,  respectively.  We  as¬ 
sessed  (computed  corresponding  values  of  0*L  for)  the 
Lyapunov  function  candidates  generated  in  step  (v)  solv¬ 
ing  Problems  1  and  2  and  further  optimized  initializing 
PENBMI  with  the  solutions  of  these  problems.  Fig.  1 
shows  0*L  and  corresponding  0*B  values  for  d{V)  =  4  and 
6.  Practically,  every  seeded  PENBMI  run  terminated 
with  the  same  0B  value  which  is  the  largest  known  (at 
least  by  us)  value  of  0  for  which  (9)  is  feasible  with  the 
prescribed  families  of  Lyapunov  functions  and  multipli¬ 
ers.  In  addition,  we  performed  10  unseeded  PENBMI 
runs  for  d(V )  =  4  and  6.  Of  these  runs  90%  and  50%, 
respectively,  terminated  successfully  (with  an  optimal 
value  of  0  equal  to  that  from  the  seeded  PENBMI  runs). 
Moreover,  unseeded  PENBMI  runs  took  longer  compu- 


Table  1 

Volume  ratios  for  ( Ei)-(£r ). 


example 

volume  ratio 

example 

volume  ratio 

(Ei) 

16.7/10.2 

(Ei) 

0.99/0.85 

(E3) 

37.2/23.5 

(E*) 

1.00/0.28 

(Eb) 

62.3  /7.3 

(Ee) 

35.0/15.3 

(Ei) 

1.44/0.70 

Fig.  1.  Histograms  of  01  (black  bars)  and  0mB  (white  bars)  sponding  references  are  contained  in  those  computed  by 

from  seeded  PENBMI  runs  for  d(V)  =  4  (left),  6  (right).  this  sequential  procedure. 


Fig.  2.  The  invariant  subsets  of  the  ROA  (dot:  d(V)  =  2, 
dash:  d(V)  =  4,  and  solid:  d{V)  —  6  (indistinguishable  from 
the  outermost  curve  for  the  limit  cycle)). 

tation  times  than  seeded  PENBMI  runs.  For  compari¬ 
son,  seeded  PENBMI  runs  took  3  —  8  and  11—24  sec¬ 
onds  for  d(V)  =  4  and  6,  respectively,  on  a  desktop  PC, 
whereas  they  took  50  —  250  and  1000  —  2500  seconds, 
respectively,  for  unseeded  PENBMI  runs.  Fig.  2  shows 
the  level  sets  of  the  Lyapunov  functions  corresponding 
to  the  value  of  0*B . 

4.2  Examples  from  the  literature 

We  present  results  obtained  using  the  method  from 
the  previous  section  for  the  systems  in  (15).  (£i)-(£3) 
are  from  (Chesi  et  al ,  2005),  (£4)  and  (£7)  are  from 
(Vannelli  and  Vidyasagar,  1985),  and  (£5)  and  (£e) 
are  from  (Hauser  and  Lai,  1992)  and  (Hachicho  and 
Tibken,  2002),  respectively.  Since  the  dynamics  in  (£1)- 
(£7)  have  no  physical  meaning  and  there  is  no  p  given, 
we  applied  SimLFG  algorithm  sequentially:  Apply 
SimLFG  algorithm  with  p(x)  =  xTx  and  Ny  =  1  for 
d(V)  =  2.  Call  the  quadratic  Lyapunov  function  ob¬ 
tained  V.  Set  p  to  V  and  apply  SimLFG  algorithm 
with  this  p  and  Ny  =  1  for  d(V)  =  4.  For  (£5)- (£7), 
we  further  applied  CWOpt  algorithm  with  Niter  =  10. 
Table  1  shows  the  ratio  of  the  volume  of  the  invari¬ 
ant  subset  of  the  ROA  obtained  using  this  procedure 
to  that  reported  in  the  corresponding  references.  Em¬ 
pirical  volumes  of  sublevel  sets  of  V  are  computed  by 
randomly  sampling  a  hypercube  containing  the  sublevel 
set.  Values  in  Table  1  axe  volumes  normalized  by  n  and 
47t/3  for  2  and  3  dimensional  problems,  respectively.  For 
(£4),  (£e),  and  (£7),  we  also  empirically  verified  that 
the  invariant,  subsets  of  the  ROA  reported  in  the  corre- 


(£1): 

(£2): 

(£3): 

(£4): 

(£5)  : 
(£«): 
(Er): 


{  ±1  =  X2,  ±2  =  —  2Xl  —  3X2  +  XjX2. 

r  x\  =  3^2? 

^  ±2  =  — 2xi  —  X2  +  X1X2  —  xf  +  X\x\  +  X%. 

(  Xi  —  X2j  X2  —  X3, 

^  X3  =  -4xi  -  3X2  -  3x3  +  x\x2  +  XjX3. 
r  Xl  =  —%2,  X2  =  — x3, 

^  x3  —  —0.915xi  +  (1  —  0.915xi)x2  —  x3. 

(  Xi  -  X2  "b  2x2X3,  X2  X3, 

^  x3  =  — 0.5xi  —  2x2  —  x3. 

{  Xl  =  -Xl  +  X2X3,  X2  =  -X2  +  X1X2,  X3  =  X3 . 

r  ii  =  -0.42xi  —  1.05x2  -  2.3Xi  —  0.5xiX2  —  x 
^  X2  =  1.98xi  +  X1X2. 

(15) 


4.3  Controlled  short  period  aircraft  dynamics 

The  closed- loop  dynamics  in  (16)  have  an  asymptotically 
stable  equilibrium  point  at  the  origin. 


E?-i  +  £<=i  ri.XiX2  +  r16x^ 

IZi=l  a2ixi  +  Y2i,j=2,5  rijxixj 

A345X 


(16) 


Here,  x  —  [xi, X2, x3, X4, X5]7,  is  the  state  vector,  xi, 
X2,  X5  are  pitch  rate,  angle  of  attack,  and  pitch  angle, 
respectively,  x3  and  X4  are  the  controller  states,  and 
A.345  G  7 Z3xb.  Before  applying  our  method,  we  performed 
excessive  simulations  and  found  a  diverging  trajectory 
whose  initial  condition  Xo  satisfies  Xq  Xo  =  16.1;  there¬ 
fore,  initialized  0s im  with  16.0.  We  applied  algorithm 
SimLFG  with  p(x)  =  xrx,  0$hrink  =  0.85,  NCOnv  = 
4000,  Ny  =  1  for  d{V)  =  2  and  4.  We  assessed  the 
Lyapunov  function  candidates  solving  Problems  A.l  and 
A. 2  and  further  optimized  using  CWOpt  algorithm  with 
Niter  =  6.  Certified  values  of  0  before  and  after  ap¬ 
plying  iterations  and  from  unseeded  PENBMI  runs  are 
shown  in  Table  2.  Unseeded  PENBMI  runs  led  to  slightly 
higher  values  of  0.  However,  this  benefit  was  at  the  ex¬ 
pense  of  high  computational  effort.  For  example,  the 
unseeded  PENBMI  rim  took  38  hours  for  d(V)  =  4 


Table  2 

Certified  values  of  0  before  and  after  applying  CWOpt  algo- 


d{V)  =  2 

d(V)  =  4 

before  iterations 

6.56 

8.99 

after  iterations 

8.56 

14.4 

PENBMI  (unseeded) 

8.60 

15.2 

Fig.  3.  A  slice  of  the  invariant  subset  of  the  ROA  (solid  line) 
and  initial  conditions  (with  =  0  and  x4  =  0)  for  diverging 
trajectories  (dots). 

whereas  our  method  took  36  minutes  (15  minutes  for  the 
SimLFG  algorithm  and  21  minutes  for  the  CWOpt  al¬ 
gorithm).  Finally,  the  dependence  that  the  starting  point 
of  CWOpt  algorithm  has  on  its  performance  is  signifi¬ 
cant.  For  example,  simply  initializing  CWOpt  algorithm 
with  V(x)  —  xT Px- h  0.001  Eti  where  PT  =  P  >-  0 
satisfies  t(P)  =  —I  yields  poor  results.  After  30  itera¬ 
tions,  the  CWOpt  iteration  converges,  but  the  resultant 
Lyapunov  function  only  certifies  Vg.s  C  Rq. 


input  and  y  is  the  output,  an  observer  L  with  polynomial 
vector  field  i  =  L(y,  z)  with  d(L)  =  3  and  a  control  law 
in  the  form  u  =  —  145.9zi  -f  12.3^2,  where  z\  and  Z2  are 
the  observer  states,  were  computed  in  (Tan,  2006).  The 
application  of  SimLFG  algorithm  with  d(V)  =  2  and  p 
from  (Tan,  2006)  and  CWOpt  algorithm  with  Niter  =  4 
lead  to  0*L  =  0.32.  We  also  applied  CWOpt  algorithm 
(initialized  with  the  quadratic  V  found  in  the  first  ap¬ 
plication)  with  d(V)  —  4  and  Nner  =  6  and  obtained 
0*L  =  0.52.  Conversely,  we  found  a  diverging  trajectory 
with  the  initial  condition  (x,  z)  satisfying  p(x,  z )  =  0.54 
proving  that  0.52  <  0*{Cl)  <  0.54. 

5  Critique  and  Conclusions 

5.1  Sampling  vs.  Simulating 

A  common  question  we  get  is  “why  simulate  to  get  the 
sample  points?  -  just  sample  some  region,  and  impose 
W(x)/(x)  <  0  there.”  There  are  a  few  answers  to  this. 
Intuitively,  even  running  a  few  simulations  gives  insight 
into  the  system  behavior.  Engineers  commonly  use  sim¬ 
ulation  to  assess  rough  measures  of  stability  robustness 
and  ROA.  Moreover,  as  converse  Lyapunov  theorems 
(Vidyasagar,  1993)  implicitly  define  a  certifying  Lya¬ 
punov  function  in  terms  of  the  flow,  it  makes  sense  to 
sample  the  flow  when  looking  for  a  Lyapunov  function 
of  a  specific  form.  Furthermore,  we  have  the  following 
observation  demonstrating  that  merely  sampling  some 
region  and  imposing  W(x)/(x)  <  0  there  carries  mis¬ 
leading  information.  Consider  the  Van  der  Pol  dynamics 
with  p(x)  =  xTx  and  let  denote  a  finite  sample  of 
Vp.  It  can  be  shown  that  the  set  of  quadratic  positive 
definite  functions  V  that  satisfy 

Si  8\{0}  C  {xelln  :  VV(x)/(x)  <  0}  (17) 


4.4  Pendubot  dynamics 

The  pendubot  is  an  underactuated  two-link  pendulum 
with  torque  action  only  on  the  first  link.  We  designed 
an  LQR  controller  to  balance  the  two-link  pendulum 
about  its  upright  position.  Third  order  polynomial  ap¬ 
proximation  of  the  closed- loop  dynamics  is  xi  =  X2, 
x2  =  782xi  -f  135x2  +  689x3  4-  90x4,  ±3  =  £4  aad 
x4  =  279x1X3  — 1425xi-257x2+273x3  — 1249x3— 171x4. 
Here,  Xi  and  x3  are  angular  positions  of  the  first  link  and 
the  second  link  (relative  to  the  first  link).  We  applied 
SimLFG  algorithm  sequentially  exactly  as  described  in 
section  4.2  and  CWOpt  algorithm  with  10  iterations  and 
obtained  0*L  =  1.69.  Conversely,  we  found  a  diverging 
trajectory  with  the  initial  condition  x  with  p(x)  =  1.95 
proving  that  1.69  <  0*{Cl)  <  1.95.  Fig.  3  shows  the 
X2  —  0  and  x4  =  0  slice  of  the  invariant  subset  of  the 
ROA  along  with  initial  conditions  (with  X2  =  0  and 
x4  =  0)  for  some  diverging  trajectories. 

4.5  Closed-loop  dynamics  with  nonlinear  observer  based 
controller 

For  the  dynamics  ±\  =  u,  x 2  =  —  £1  +  x?/6  —  u  and 
y  =  x 2,  where  xi  and  X2  are  the  states,  u  is  the  control 


is  nonempty.  In  fact,  for  V(x)  =  0.32xf  —  0.25xiX2  + 
0.31x2,  (17)  is  satisfied  (actually  for  all  x  €  Vi  s, 
W (x)f(x)  <  —  f3(x)).  This  naively  suggests  to  draw 
samples  from  the  set  of  quadratic  positive  definite 
functions  satisfying  (17)  in  order  to  try  to  prove  that 
V\.s  C  Ho.  However,  simulations  reveal  a  contradicting 
fact:  Using  trajectories  with  initial  conditions  in  Si.s  for 
d(V)  =  2,  i.e.,  with  tp(x)  =  [xf,  £1X2,  £i]T?  constraints 
(6),  (11)  (with  7  =  1),  and  (12)  turn  out  to  be  infeasi¬ 
ble.  This  verifies  that  no  quadratic  Lyapunov  function 
can  prove  V\.g  C  Ro  using  conditions  (6)-(8),  with  the 
additional  constraint  that  V{x)  <  —  10~6xrx  on  all 
trajectories  starting  in  Vi  s-  Recall  though,  that  using 
quartic  Lyapunov  functions  we  know  0*  (V^S)  >  2.14. 
By  these  observations,  we  have  the  following  series  of  in¬ 
clusions  for  the  subsets  of  the  positive  definite  quadratic 
polynomials 

{V  :  V  certifies  V$  C  Ro  using  (6)  -  (8)} 

C  {V  :  W(cs(r))/(cs(r))  <  0  Vr,  Vs  €  S*} 

C  {V  :  W(s)/(s)  <  0  V  s  G  S*}, 


where  cs  denotes  the  trajectory  with  the  initial  condi¬ 
tion  s  €  S Therefore,  merely  sampling  instead  of  using 
simulations  leads  to  a  larger  outer  set  from  which  the 
samples  for  V  are  taken  in  step  (v)  of  SimLFG  algo¬ 
rithm  and  it  is  less  likely  to  find  a  function  that  certifies 
that  C  i2o- 

5.2  Conclusions 

We  proposed  a  method  for  computing  invariant  subsets 
of  the  region-of- attraction  for  asymptotically  stable  equi¬ 
librium  points  of  dynamical  systems  with  polynomial 
vector  fields.  We  used  polynomial  Lyapunov  functions  as 
local  stability  certificates  whose  certain  sublevel  sets  are 
invariant  subsets  of  the  region-of-attraction.  Similar  to 
many  local  analysis  problems,  this  is  a  nonconvex  prob¬ 
lem.  Furthermore,  its  sum-of-squares  relaxation  leads  to 
a  bilinear  optimization  problem.  We  developed  a  method 
utilizing  information  from  simulations  for  easily  gener¬ 
ating  Lyapunov  function  candidates.  For  a  given  Lya¬ 
punov  function  candidate,  checking  its  feasibility  and 
assessing  the  size  of  the  associated  invariant  subset  are 
affine  sum-of-squares  optimization  problems.  Solutions 
to  these  problems  provide  invariant  subsets  of  the  region- 
of-attraction  directly  and/or  they  can  further  be  used 
as  seeds  for  local  bilinear  search  schemes  or  iterative 
coordinate- wise  linear  search  schemes  for  improved  per¬ 
formance  of  these  schemes.  We  reported  promising  re¬ 
sults  in  all  these  directions. 
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A  Appendix 

Lemma  2  Given  goi9u"  ,9m  €  R[x],  if  there,  exist 
$i,  ■  •  •  ,  $m  G  £[x]  such  that  go  —  Yl’iLi  si9i  €  E[x],  then 
{x  €  Un  :  g\{x),.. ., gm{x )  >  0}  C  {x  €  7Zn  :  g0(x)  >  0}  . 

Lemma  3  Given  <?i  >  02  €  R[x]  such  that  g§  is  positive 
definite  and  go{0)  =  0,  if  there  exist  51,32  €  £[x]  such 
that  g\Sx  4-  92S2  ~9o  €  E[x],  then  {x  G  Tln  :  <71  (x)  < 
0}\{0}  C  {x  G  Tin  :  92{x)>  0}.  « 

Problems  1  and  2  in  section  3  compute  lower  bounds 
on  the  largest  value  of  7  and  (3  such  that,  for  given  V 
and  p,  fivvA{0}  C  {x  G  7ln  :  W(x)f(x)  <  0}  and 
V(3  C  We  propose  alternative  formulations,  that  do 
not  require  line  search,  to  compute  similar  lower  bounds. 
Labeled  7*  and  /?*,  these  are  generally  different  than  7 *L 
and  /?£.  For  h,g  G  Ri]  and  a  positive  integer  d,  define 


fi°(h,  g )  :=  infz/o  Hx)  such  that  g(x)  =  0,  and 


fi*(h,g,d)  :=  sup  p  subject  to 

/x>0,r€R[x] 

(*  -  m)  (a:?'1  +  •  •  •  +  -  gr  e  E[x]. 

Note  that  jj,*(h,g,d)  <  p°(h,g). 

Lemma  4  Let  g,  h  :  7£n  — >  TZ  be  continuous ,  /i 
6e  positive  definite,  g( 0)  =  0,  and  g(x)  <  0  /or  all 
nonzero  x  €  O,  a  neighborhood  of  the  origin.  De¬ 
fine  7°  :=  p°(h,g).  Then,  the  connected  component  of 
{x  £  7 Zn  :  h(x)  <  7°}  containing  the  origin  is  a  subset 
of  {x  £  7Zn  :  g(x)  <  0}  U  {0}.  <i 

Proof:  Suppose  not  and  let  z  ^  0  be  in  the  con¬ 
nected  component  of  {z  £  1 Zn  :  h(x)  <  70}  containing 
the  origin  but  g(x)  >  0.  Then,  there  exists  a  contin¬ 
uous  function  d  :  [0, 1]  — ►  7Zn  such  that  t?(0)  =  0, 
t?(l)  =  x,  and  h(d(t))  <  70  for  all  t  £  [0,1].  Since 
<?(0)  =  0  and  g(x)  <  0  for  all  nonzero  x  £  O,  there 
exists  0  <  6  <  1  such  that  <?(t?(e))  <  0.  Since  z  is 
not  in  {z  £  7 Zn  :  g(x)  <  0},  y(t?(l))  >  0.  Since  g  and 
d  are  continuous,  there  exists  t*  £  (0, 1]  such  that 
=  0,  which  implies  h(d(tm))  >  70.  This  contra¬ 
diction  leads  to  z  £  {z  £  7 Zn  :  g(x)  <  0}.  □ 

Corollary  5  Let  V  £  R[z]  be  a  positive  definite  Cl  func¬ 
tion  and  satisfy  (12)  andV( 0)  =  0.  Then,  for  all  7  such 
that  0  <  7  <  p0{V,W f),  the  connected  component  of 
fly, 7  containing  the  origin  is  an  invariant  subset  of  the 
ROA.  < 

Proof:  Since  the  quadratic  part  of  V  is  a  Lyapunov 
function  for  the  linearized  system,  there  exists  a  neigh¬ 
borhood  O  of  the  origin  such  that  W(z)/(z)  <  0 
for  all  nonzero  z  £  O.  By  Lemma  (4),  the  connected 
component  of  fly, 7  containing  the  origin,  a  subset  of 
the  connected  component  of  {z  £  7Zn  :  V(x)  < 
p°(V,  VVr/)}  containing  the  origin,  is  contained  in 
{z  £  7Zn  :  W(x)f(x)  <  0}U{0}.  Corollary  (5)  follows 
from  regular  Lyapunov  arguments  (Vidyasagar,  1993). 
□ 

Corollary  6  For  some  positive  integer  d\,  define  7*  := 
p*(V,  W/,  di).  Then ,  if  7  <  7*  for  some  positive  integer 
d\,  then  the  connected  component  of  tty, 7  containing  the 
origin  is  an  invariant  subset  of  the  ROA .  < 

Corollary  7  Let  0  <  7  <  7*,  d%  be  a  positive  integer, 
V,p  £  R[z]  be  positive  definite  andp  be  convex.  Define 
Pi  •=  V'-7,d2).  Then  for  any  0  <  01,  Vp  C 

and  VpCRo-  < 


Stability  Region  Analysis  using  polynomial  and 
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Abstract 

We  propose  using  (bilinear)  sum-of-squares  programming  for  obtaining  inner  bounds  of  regions- 
of-attraction  for  dynamical  systems  with  polynomial  vector  fields.  We  search  for  polynomial  as  well  as 
composite  Lyapunov  functions,  comprised  of  pointwise  maximums  of  polynomial  functions.  Results  for 
several  examples  from  the  literature  are  presented  using  the  proposed  methods  and  the  PENBM1  solver. 


I.  Introduction 

Finding  the  stability  region  or  region-of-attraction  (ROA)  of  a  nonlinear  system  is  a  topic 
of  significant  importance  and  has  been  studied  extensively,  for  example  in  [1-12].  It  also  has 
practical  applications,  such  as  determining  the  operating  envelope  of  aircraft  and  power  systems. 

Most  computational  methods  aim  to  compute  an  inner  bound  on  the  region-of-attraction, 
namely  a  set  that  contains  the  equilibium  point,  and  is  contained  in  the  region-of-attraction. 
The  methods  above  can  roughly  be  split  into  Lyapunov  and  non-Lyapunov  methods.  Lyapunov 
methods  (the  focus  of  this  paper)  are  based  on  local  stability  theorems  and  search  for  functions 
satisfying  conditions  which  quantitatively  prove  local  stability.  Nonlinear  programming  is  used  in 
[1]  to  optimize  (by  choice  of  positive  definite  matrix)  the  volume  of  an  ellipsoid  contained  in  the 
region-of-attraction.  Rational  Lyapunov  functions  that  approach  oo  on  the  boundary  of  the  region- 
of-attraction  are  constructed  iteratively  in  [2],  motivated  from  Zubov’s  work.  Computational  con¬ 
siderations  limit  the  degree  of  the  rational  function,  and  inner  estimates  to  the  ROA  are  obtained. 
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Easy  to  compute  estimates  are  considered  in  [3],  which  restricts  the  Lyapunov  function  search  to 
low-dimensional  manifold  of  quadratic  Lyapunov  functions,  obtaining  analytical  simplifications. 
Following  [1],  but  employing  semidefinite  programming  techniques,  [4]  aims  to  maximize  the 
volume  of  an  ellipsoid  whose  containment  in  the  region-of-attraction  can  be  ascertained  with 
sum-o f-squares  (SOS)  decompositions.  Attention  is  restricted  to  odd,  polynomial  vector  fields, 
and  SOS  optimization  is  combined  with  general  nonlinear  programming.  A  sequence  of  functions, 
called  nested  Lyapunov  functions,  are  introduced  in  [5]  to  derive  stability  region  and  rate- 
of-convergence  estimates.  Both  [6]  and  [7]  solve  a  sequence  of  linear  semidefinite  programs, 
iteratively  searching  over  Lyapunov  function  candidates  and  sum-of-squares  multipliers.  The 
“coordinatewise”  ascent  method  is  generally  effective,  though  no  convergence  result  holds.  By 
contrast,  the  formulation  here  is  more  direct,  but  yields  a  single  bilinear  (nonconvex)  SOS 
program.  Closely  related  to  Lyapunov  methods  are  viability  methods,  which  effectively  integrate 
an  invariant  set  backwards  in  time,  obtaining  increasing  estimates  for  the  region-of-attraction. 
Both  [8]  and  [9]  use  discretization  (in  time)  to  flow  invariant  sets  backwards  along  the  flow 
of  the  vector  field,  obtaining  larger  and  larger  estimates  for  the  region-of-attraction.  In  [8],  the 
invariant  sets  are  restricted  to  be  sublevel  sets  of  polynomials,  and  the  discretized  backwards  flow 
is  approximated  with  semidefinite  programming.  The  approach  of  [9]  also  requires  discretization 
in  space  and  suffers  from  exponential  growth  in  state  dimension.  Generally,  the  method  is  exact, 
but  computation  may  require  exponential  growth  in  state  dimension.  Depending  on  the  user’s 
point  of  view,  problems  of  modest  (between  4  and  8)  state  dimension  are  intractable.  Non- 
Lyapunov  methods,  [10]  and  [11]  focus  on  topological  properties  of  regions  of  attraction.  A 
survey  of  results,  as  well  as  an  extensive  set  of  examples  and  new  ideas,  is  presented  in  [12]. 

In  this  paper,  we  present  a  method  of  using  sum-of-squares  (SOS)  programming  to  search 
for  polynomial  Lyapunov  functions  that  enlarge  an  inner  estimate  of  the  region-of-attraction  of 
nonlinear  systems  with  polynomial  vector  fields.  SOS  programming,  coupled  with  polynomial 
Lyapunov  functions  has  roots  that  can  be  traced  back  at  least  to  Bose  and  Li  [13]  and  Brockett, 
[14]  and  the  power  transform  of  Barkin  et.al  [15],  which  was  used  in  [16]  to  find  non-quadratic 
Lyapunov  functions  for  uncertain  linear  systems.  Recent  theoretical  work,  [17],  [7]  and  [18], 
continues  to  further  the  role  of  this  approach.  An  impediment  to  using  high  degree  Lyapunov 
functions  is  the  extremely  rapid  increase  in  the  number  of  optimization  decision  variables  as  the 
state  dimension  and  the  degree  of  the  Lyapunov  function  (and  the  vector  field)  increase.  Here, 


we  propose  using  pointwise  maximums  of  polynomial  functions  to  obtain  rich  functional  forms 
while  keeping  the  number  of  optimization  decision  variables  relatively  low.  Pointwise  maximum 
and  other  composite  Lyapunov  functions  have  been  used  in  many  instances,  [19],  [20],  [21], 
including  stability  and  performance  analysis  of  constrained  systems  and  robustness  analysis  of 
uncertain  systems,  where  affine  and  polynomial  parameter-dependent  Lyapunov  functions  are  also 
used,  [22],  [23].  The  notation  is  generally  standard,  with  IZ^  denoting  the  set  of  polynomials 
with  real  coefficients  in  n  variables  and  C  Tin  denoting  the  subset  of  SOS  polynomials. 

II.  Estimating  a  Region  of  Attraction 
Consider  an  autonomous  dynamical  system  of  the  form 

x{t)  =  f(x(t ))  (1) 

where  x{t)  G  Rn  and  /  is  an  n-vector  of  elements  of  TZn  with  /( 0)  =  0.  The  following  lemma 
on  invariant  subsets  of  the  region-of-attraction  is  a  modification  of  ideas  from  [24,  pg.  167]  and 
[25,  pg.  122]: 

Lemma  I:  If  there  exist  continuously  differentiable  functions  {K}^=1  :  Rn  — *•  E  such  that 


V(x )  :=  max  VAx)  is  positive  definite,  (2) 

\<i<q 

ft  :=  {x  G  |  V(x)  <  1}  is  bounded,  (3) 

Li  :=  (i  €  ln  |  max  Vj(x)  <  Vfx)  <  1},  i  =  1, . . . ,  q  (4) 

i  <j<q 

Li  \  {0}  C  {x  G  Rn  I  ^/(x)  <  0},  i  =  1, . . . ,  q,  (5) 


then  for  all  x(0)  G  ft,  the  solution  of  (1)  exists,  satisfies  x(t)  G  ft,  and  limt^oox(t)  =  0.  As 
such,  ft  is  invariant,  and  a  subset  of  the  region-of-attraction  for  (1). 

Proof:  The  proof  is  written  for  q  =  2.  The  extension  to  q  =  1  or  q  >  2  is  straightforward. 
Since  L^  U  L2  =  ft,  condition  (5)  insures  that  if  x(0)  G  ft,  V(x(t))  <  V(x(0))  <  1  while  the 
solution  exists.  Solutions  starting  in  ft  remain  in  ft  while  the  solution  exists.  Since  ft  is  compact, 
the  system  (1)  has  an  unique  solution  defined  for  allt  >  0  whenever  x(0)  G  ft. 

Take  e  >  0.  Define  Se  (x  G  Rn  1 1  <  V(x)  <  1},  so  St  C  (Lj  U  Lf)  \  {0}.  Note 
that  for  each  i,  (Sc  D  L{)  C  Li  \  {0}  C  {x  G  En  |  ^/(x)  <  0},  so  on  the  compact  set 
S'e  D  Lu  3  rit£,  such  that  ^/(x)  <  — rii£  <  0.  Consequently,  if  x(t)  G  Se  D  L\  on  [ tA,tB ], 


then  V(x(tB ))  <  -ru(tB  -  tA)  +  V(x(tA)).  Similarly,  if  x(t)  €  S£  n  L2  on  [tA,tB],  then 
V(x(tB ))  <  -r2,e(tB  -  tA)  +  V(x(tA)).  Therefore,  if  x(t)  €  S£  D  (Li  U  L2)  on  [tA,tB],  then 
V(x(tB))  <  -re(tB -tA)  +  V(x(tA)),  where  r£  =  min(rij£,r2,£).  Since  r£  >  0,  this  implies  that 
3t*  such  that  I/(a:(f))  <  e  for  all  t  >  t*,  i.e.  x(t)  €  T£  :=  {i  G  1”  |  V(x)  <  e}  for  all  t  >  t*. 
This  shows  that  if  x(0)  €  (2,  V(x(t))  — ►  0  as  t  — *  oo. 

Finally,  let  e  >  0.  Define  Qe  :=  {x  €  Rn  |  ||x||  >  e,  V^a;)  <  1}.  J2£  is  compact,  with  0  ^  ile. 
Since  V  is  continuous  and  positive  definite,  37  such  that  V(x)  >  7  >  0  on  ile.  We  have  already 
established  that  — *  0  as  t  — *  00,  so  3t  such  that  for  all  t  >  t,  V(a:(t))  <  7  and  hence 

x(t)  Q£,  which  means  ||a:(t)||  <  e.  So  x(t)  — *  0  as  t  — ►  00.  ■ 

Remarks:  Standard  modifications  to  the  hypothesis  of  Lemma  1  can  yield  global  stability 
conditions  as  well.  However,  neither  formulation  can  yield  exact  results  for  systems  whose 
region-of-attraction  is  unbounded,  but  not  all  of  Rn  (since  in  Lemma  1,  il  must  be  is  bounded). 
See  section  III-C  for  further  details.  The  constraints  in  equations  (2)-(5)  are  not  convex  constraints 
on  V,  as  illustrated  by  a  1 -dimensional  example,  [26].  Take  f(x)  =  —x,  q  =  1  and  = 

16x2  -  19.95x3  +  6.4a:4  and  V?(x)  =  0.1a:2.  Then  Vf  and  V?  satisfy  (2)-(5),  but  0.58V}1  +  0.42V,6 
does  not. 

In  order  to  enlarge  £2  (by  choice  of  V),  we  define  a  variable  sized  region  P$  :=  (a:  G 
Rn  |  p(x)  <  0},  and  maximize  0  while  imposing  the  constraint  Pp  C  F2.  Here,  p(x)  is  a  fixed, 
positive  definite  polynomial,  chosen  to  reflect  the  relative  importance  of  the  states.  Applying 


Lemma  1,  the  problem  is  posed  as  an  optimization: 

max  0  s.t.  V^(0)  =  0 
pm 

V(x)  :=  max  Vi(x)  is  positive  definite,  (6) 

1  <i<q 

f2  :=  {x  €  Rn  |  V”(a:)  <  1}  is  bounded,  (7) 

Pp  c  n  (8) 

{x  €  Rn  |  max  Vs(x)  <  V^x)  <  1}  \  {0}  C  {x  e  Rn  |  f  }f(x)  <  0}  (9) 


where  (9)  holds  for  i  —  1  Let  l\(x)  be  a  fixed,  positive  definite  polynomial.  For  each  Vt, 

if  we  require  Vi—li  €  £n  for  i  =  1, . . . ,  q,  then  both  (6)  and  (7)  are  satisfied.  Clearly,  (8)  holds 


if  and  only  if 


Q 

{x  e  Rn  |  p(x)  <(3}  C  f]{x  G  Rn  |  Vi(x)  <  1},  (10) 

i= 1 

Introducing  another  fixed,  positive  definite  polynomial,  l2(x),  we  can  apply  Lemmas  2  and  3 
(see  appendix)  to  obtain  sufficient  conditions  which  ensure  constraints  (9)  and  (10)  hold.  Written 
as  an  optimization,  the  problem  is 

max [3  over  3  G  R,  G  Vt(0)  0)  $2ti  ^3ti  ® oij  G  •  •  •  Q 

such  that 

Vx  —  h  £n,  (11) 

—  {(P  ~  P)su  +  (K  —  1)  j  £  £n,  (12) 

Q 

~  [(1  -  V{)s2i  +  gf/S3i  +  i2]  -  e  En-  (13) 

j= i 
}*i 

All  constraints  are  sum-of-square  constraints,  however  (even  for  q  =  1)  products  of  decision 
variables  are  present.  Therefore,  the  optimization  cannot  be  translated  into  a  linear  semidefinite 
program,  but  is  converted  to  a  bilinear  semidefinite  program.  Two  of  the  conditions  require 
positivity  (beyond  nonnegativity),  and  the  fixed  positive-definite  polynomials,  l\  and  l2  are 
introduced  as  offsets  to  enforce  this.  Next  we  present  results  from  several  small  problems.  We 
have  chosen  to  rely  on  the  PENBMI  solver  [27],  a  local  bilinear  matrix  inequality  solver  from 
PENOPT  to  attack  our  problems.  This  uses  a  penalty  method.  Alternate  approaches  to  BMIs, 
such  as  linearization  and  homotopy,  [28]  and  interior  point  methods,  [29,  Chap  7],  may  yield 
improved  results  and/or  superior  computational  efficiency.  Resolving  these  questions  is  left  for 
further  research. 


III.  Examples 

All  of  the  systems  considered  are  locally  exponentially  stable.  The  notation  ny  denotes  the 
degree  of  V,  specifically  each  Vi  includes  all  monomials  of  degree  2  up  through  ny.  In  all 
examples,  p  is  quadratic,  and  the  degree  of  su  is  chosen  so  that  the  degree  of  the  polynomial 
in  equation  (12)  is  equal  to  ny.  The  integer  nA  denotes  degree  of  the  polynomial  in  (13).  Once 
ny  is  chosen,  and  the  vector  field  /  is  fixed,  nA  limits  the  degrees  of  the  multipliers  s2i,  s2l 
and  SQij  through  simple  degree  counting.  In  each  case,  the  positive  definite  polynomials  l\  and 


k  are  of  the  form  lk(x)  =  Yh=i  ek<ix j.  For  the  purposes  of  computation,  the  efc,t  are  treated  as 
additional  decision  variables,  and  constrained  to  satisfy  ekti  >  10~7. 

A.  Example  1  -  Van  der  Pol  equations 

The  system  is  ii  =  — x2,x2  =  X\  +  {x\  —  l)x2.  It  has  an  unstable  limit  cycle  and  a  stable 
equilibrium  point  at  the  origin.  Finding  its  region-of-attraction  has  been  studied  extensively,  for 
example,  in  [1],  [12],  [11].  The  region-of-attraction  for  this  system  is  the  region  enclosed  by  its 
limit  cycle,  which  is  easily  visualized  from  the  numerical  solution  of  the  ODE.  However,  our 
goal  is  to  use  the  bilinear  SOS  formulation.  For  this  example,  p  is  chosen  to  be  xTRx,  for  two 
different  €  R2x2, 


0.38 

-0.14 

,  and  f?2  := 

0.28 

0.14 

-0.14 

0.28 

0.14 

0.38 

The  results  using  shape  factor  defined  by  R\  using  the  pointwise  maximum  of  two  fixed 
degree  polynomials  are  listed  in  Table  I.  Fig.  1  shows  the  limit  cycle  and  the  level  sets  of 
the  certifying  Lyapunov  functions1.  The  level  set  of  pointwise  maximum  of  two  6th  degree 
polynomial  functions  includes  nearly  the  entire  actual  region-of-attraction.  The  dashed  line  is 
the  level  set  of  p  (for  ny  =  6),  which  clearly  shows  that  our  p  has  been  preselected  to  “align” 
closely  with  the  actual  region-of-attraction.  Of  course,  this  would  be  impossible  to  do  in  general, 
and  we  discuss  the  implications  of  this  later  in  this  section.  Our  results  compare  favorably  with 
[11]  as  well  as  the  degree  6  solution  from  [7],  and  the  final  (40th)  iterate  from  degree  6  solutions 
of  [8],  all  of  which  are  shown  in  Figure  2.  Clearly,  the  solution  of  [8]  is  a  very  high  quality 
estimate  of  the  true  ROA.  Parametrizing  the  boundaries  using  polar  coordinates  reveals  that  as 
a  function  of  angle,  the  radius  of  [8]  exceeds  our  ny  =  6  radius  on  52.6%  of  [0  2tt];  is  0.22% 
larger,  on  average,  than  our  ny  =  6  radius;  exceeds  our  ny  =  6  radius  by  as  much  as  1 .4%  in 
some  directions;  is  smaller  than  our  ny  =  6  radius  by  as  much  as  0.8%  in  other  directions.  We 
conclude  that  the  result  in  [8]  is  very  similar,  though  slightly  superior  to  our  result. 

It  is  interesting  to  observe  how  the  V  functions  interact  in,  for  example,  the  6th  degree  case 
Figure  3  shows  the  level  sets  {x\Vi(x)  <  1}.  For  VJ,  there  are  3  connected  components,  one 


lThe  certifying  Lyapunov  functions  and  SOS  multipliers  for  all  examples  in  this  paper  are  available  at 
http : // jagger .me .berkeley . edu/~pack/certif icates 


“large”  component  centered  at  the  origin  (whose  boundary  is  essentially  the  limit  cycle),  and 
2  “islands”  in  the  2nd  and  4th  quadrants.  For  V2,  the  level  set  is  one  connected  component 
centered  at  the  origin,  visually  the  same  as  the  large  component  of  Vi .  Label  the  two  islands  as 
Ii  and  /2,  and  the  intersection  of  the  two  (nearly  identical)  centered  components  as  il. 

Inside  f  and  I2,  V\  it  0,  but  V2  >  1  >  Vj,  so  I\  and  I2  are  excluded  in  the  set  ft.  Moreover, 
on  Q,  Vi  <  0  where  Vi  >  V2,  and  V2  <  0  where  V2  >  Vi,  proving  that  Cl  is  a  region-of-attraction. 
Since  {x  |  V2(x)  <  1}  «  il,  it  is  tempting  to  assume  that  V2  alone  can  prove  the  stability  claim. 
However,  many  points  have  V2  >  0  (the  shaded  region  in  il). 

In  this  example,  using  pointwise  maximum  of  three  polynomials  does  not  offer  additional 
benefits  (row  1  and  4  of  Table  I).  Better  results  are  obtained  (row  5)  by  increasing  the  degree 
of  the  {sj},  but  this  increases  the  number  of  decision  variables,  so  the  computational  benefit  is 
effectively  erased. 

Finally,  optimizing  with  the  shape  factor  defined  by  i?2  yields  almost  identical  results  (in 
terms  of  17).  Fig.  4  illustrates  the  analogous  level  sets  of  V,  and  also  shows  a  level  set  for  this 
p.  Clearly,  the  level  sets  for  this  shape  factor  are  not  aligned  with  the  actual  region-of-attraction, 
nevertheless,  the  optimization  performs  quite  well. 

B.  6  examples  from  [4] 

Reference  [4]  aims  to  maximize  the  volume  of  an  inner  ellipsoidal  estimate  of  the  region  of 
attraction,  presenting  results  from  6  examples.  The  volume  reported  in  [4]  is  normalized:  in 
R2  it  is  2-dimensional  area  divided  by  7 r,  while  in  R3  it  is  3 -dimensional  volume  divided  by 
As  an  exercise,  we  solve  the  same  problems  here.  The  results  are  summarized  in  Table  II. 
Maximizing  volume  is  not  directly  compatible  with  our  scalar  objective  involving  the  function 
p  (whose  level  sets  may  or  may  not  be  ellipsoidal).  We  began  with  a  simple  approach:  using 
a  spherical  shape  factor,  p(x)  :=  xTx,  solve  the  optimization  problem  and  then  compute  the 
volume  of  the  level  set  {x  :  V(x)  <  1}  (easily  computed  for  a  quadratic  V,  and  estimated  with 
Monte  Carlo  integration  for  high  degree  and  pointwise-max  V’s).  Problems  SI,  S2,  S3  and  S4 
are  successfully  addressed  using  this  approach.  Note  the  improvement  for  S2  when  the  degree 
of  the  multipliers  is  increased  (via  nf)  even  though  ny  is  held  constant.  Problem  S5  required 
an  alteration,  referred  to  as  bootstrap,  to  obtain  large  volumes.  In  this  calculation,  the  initial 
optimization  was  as  above,  with  a  spherical  p,  using  quadratic  Lyapunov  function  candidates. 


Subsequent  optimizations,  with  richer  Lyapunov  function  candidates  used,  for  p,  the  obtained 
quadratic  Lyapunov  function  (as  opposed  to  xTx).  Problem  S6  is  more  challenging  and  the 
methods  we  present  here  do  not  obtain  volumes  as  large  as  those  reported  in  [4],  The  S6  table 
entry  involving  quartic  functions  is  empty,  as  PENBMI  exhibited  unreliable  behavior  on  this 
problem,  exposing  some  genuine  deficiences  in  our  overall  approach. 

C.  Unbounded  Region-of-Attraction 

Consider  X\  =  x2,x2  =  -(1  -  x\)x\  -  x2  from  [30].  The  region  of  attraction  to  the  stable 
equilibrium  at  x  =  0  is  unbounded,  but  not  all  of  M2.  Exact  methods,  such  as  those  in  [9],  may 
obtain  the  correct  answer  in  this  problem.  By  contrast,  the  formulation  in  equations  (11)-(13) 
cannot,  since  il  is  necessarily  compact.  Using  a  simple  p(x)  :=  x\  +  x\  shape  factor,  we  obtain 
nearly  identical  results  for  quadratic  and  pointwise-max  quadratic  Lyapunov  functions,  yielding 
P  such  that  Pp  nearly  touches  the  stability  boundary,  and  the  bounded  level  set  {x  :  V(x)  <  1}  is 
ellipsoidal,  roughly  aligned  with  the  true  region-of-attraction.  Using  the  bootstrap,  with  nv  =  6 
yields  significant  improvement.  The  two  level  sets  are  shown  in  the  left  panel  of  figure  5,  along 
with  some  trajectories  of  the  system. 

D.  An  example  from  reference  [2] 

Another  2-state  example  with  polynomial  vector  field  comes  from  example  4  in  [2].  The 
dynamics  are  x\  =  — 0.42^!  —  1.05x2  —  2.3x2  —  0.5xiX2  —  x\\x2  =  1.98xi  +  x\x2.  The  inner 
estimate  from  [2]  along  with  our  estimate  using  quadratic,  quartic,  pointwise-max  quartic,  and 
degree  6  functions  are  shown  in  the  right  panel  of  Figure  5.  Pointwise-max  ( q  =  2)  degree  6 
solutions  yielded  no  appreciable  improvement  over  the  q  =  1  case,  and  are  not  shown. 

IV.  Benchmark  Study 

There  are  several  drawbacks  to  our  approach,  most  notably  searching  over  the  non-convex 
decision  variable  space.  Given  this  deficiency,  it  is  useful  to  investigate  how  equations  (11)- 
(13),  coupled  with  the  PENBMI  solver  perform  on  an  “easy”  nonlinear  problem,  with  respect 
to  “arbitrary”  data  and  increasing  problem  size.  Let  x  =  —lx  +  (xTBx)x  where  x(t)  £  Rn,  and 
B  £  Rnxn,  B  >-  0.  For  this  system,  inspired  by  Example  5  of  [1],  the  set  {x  £  Rn  |  xTBx  <  1} 
is  the  exact  region-of-attraction  for  the  x  =  0  equilibrium  point  (use  V (x)  :=  xTBx  to  prove 


this).  Let  Pp  :=  {x  E  Rn  |  xTRx  <  /3},  R  E  Mnxn,  R  >-  0.  The  supremum  value  for  (3  such 
that  Pp  C  {x  E  P.n  |  xTBx  <  1}  is  (3  =  [\max{R~1BR~^)]~1.  Equations  (11  )-(l 3)  can  yield 
this  answer,  specifically,  take  q  =  1  and  for  any  7  >  1,  choose  r  such  that  1  <  r  <  7.  Then 
for  large  enough  a  (depending  on  fixed  choice  of  quadratic  l2)  the  choices  V(x)  :=  7 xTBx, 
s2  :=  2 arxTBx  and  s3  :=  a  satisfy  (13),  prove  that  {x\xTBx  <  1}  is  in  the  region-of- 
attraction.  Hence,  this  class  of  problems  provides  a  test  for  any  specified  BMI  solver  to  actually 
discover  the  known-to-exist  solution.  For  each  n,  100  trials  are  performed.  Each  trial  consists  of 
a  random  choice  of  positive  definite  B  and  R,  each  with  eigenvalues  exp(2rt)  where  each  r,  is 
picked  from  a  normal  distribution  with  zero  mean  and  unit  variance,  and  random,  orthonormal 
eigenvectors.  For  each  trial,  we  run  the  PENBMI  optimizer  3  times  (initial  point  randomly  chosen 
each  run).  Table  III  shows  the  results  of  the  test. 

A  run  is  classified  successful  if  the  solver  returns  the  message  “No  problems  detected”,  and 
classified  failure  otherwise.  Except  for  the  case  of  n  =  6,  there  are  no  trials  that  fail  for  all  3  runs 
(for  n  =  6,  one  trial  did  fail  in  all  3  runs,  and  note  that  this  single  instance,  3-trial  failure  is  not 
taken  into  account  in  the  table  entries  described  below).  Among  the  successful  runs,  the  quality 
of  the  answer  is  assessed  by  the  nearness  of  (3x  Amax  to  1.  The  worst  case  (smallest)  value  among 
the  (296-300)  successful  runs  is  given.  The  next  column  shows  the  worst  case  (3  x  Amax  over 
100  trials,  exploiting  the  3  repeated  attempts  and  the  randomized  initial  starting  point  chosen  by 
PENBMI.  The  entries  are  «  1,  which  indicates  that  repeated  runs  of  the  same  problem  eventually 
lead  to  the  optimal  solution  for  this  example.  For  this  limited  benchmark  example,  although  our 
problem  formulation  is  bilinear  in  the  decision  polynomials  and  the  bilinear  solver,  PENBMI, 
is  a  local  solver,  the  results  obtained  are  encouraging. 

V.  Conclusions 

In  this  paper,  we  presented  techniques  using  sum-of-squares  programming  for  finding  provable 
regions-of-attraction  for  nonlinear  systems  with  polynomial  vector  fields.  Several  small  examples 
are  presented.  For  systems  with  cubic  vector  fields,  analyzing  local  stability  using  Lyapunov 
functions  which  are  the  pointwise-max  of  quadratic  and  quartic  functions  appears  to  be  a  useful, 
and  modestly  tractable  extension  to  simply  using  polynomial  Lyapunov  functions. 
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VII.  Appendix 

A  monomial  ma  in  n  variables  is  a  function  defined  as  ma{x )  =  xa  :=  f°r 

q  €  Z".  The  degree  of  a  monomial  is  defined,  deg ma  ^”=1  a,.  A  polynomial  /  in  n 
variables  is  a  finite  linear  combination  of  monomials,  with  ca  G  R: 

/:=^cQmQ  =  ^cQx“. 

a  a 

Define  IZn  to  be  the  set  of  all  polynomials  in  n  variables.  The  degree  of  /  is  defined  as  deg  /  := 
maxQdegmQ  (provided  the  associated  cQ  is  non-zero).  Additionally  define  En  to  be  the  set  of 
sum-of-squares  (SOS)  polynomials  in  n  variables. 

t 

P=YJ?  ,  fi  eTZn,  **  1, 

i=l 

Obviously  if  p  €  En,  then  p(x)  >  0  Vrr  €  Rn.  A  polynomial,  p  €  En  iff30^(5€Krxr  such 
that  p(x)  =  zT(x)Qz(x),  with  z(x)  a  vector  of  suitable  monomials.  The  set  of  Q  that  satisfies 
p(x)  =  zT(x)Qz(x)  is  an  affine  subspace,  so  that  semidefinite  programming  plays  the  key  role 
in  deciding  if  a  given  polynomial  is  SOS.  The  lemmas  below  are  elementary  extensions  of  the 
<S-procedure,  [32],  and  very  limited  special  cases  of  the  Positivstellensatz,  [33,  Theorem  4.2.2] 
In  both  cases,  the  SOS  polynomials  {sjt}^!  are  often  called  the  “SOS  multipliers.” 

Lemma  2:  Given  Pi,P2  €  Tin,  and  positive  definite  h  €  7 Zn,  with  h( 0)  =  0.  If  Si,  s2  €  En 
satisfy  piSi  +P2S2  —  h  €  E  then  {x  :  P\{x)  <  0}  \{0}  C  {x  :  pi{x)  >  0}. 


Sn  :=  {  P  €  Tin 


Lemma  3:  Given  {p,}£o  €  7 In.  If  there  exist  {s*}^  €  En  such  that  pQ  -  J2?=i  siPi  €  S„, 
then  f)T=l{x  €  En  | p,(x)  >  0}  C  {x  €  Kn  | po(x)  >  0}. 

SOSTOOLS,  [34],  [35],  GloptiPoly,  [36],  and  YALMIP,  [31]  automate  the  translation  from 
SOS  programs  to  semidefinite  programs,  converting  to  solver-specific,  e.g.,  SeDuMi  [37]  or 
SDPT3  [38],  syntax.  YALMIP  also  handles  bilinear  decision  polynomials,  using  PENBMI  [27]. 

Despite  these  software  tools,  and  even  ignoring  the  nonconvexity  of  our  formulation,  there 
are  significant  dimensionality  problems  as  well:  [39,  Table  6.1]  illustrates  the  unpleasant  growth 
in  the  number  of  decision  variables  with  n  and  the  polynomial  degree. 
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table  i 

Pointwise  Max  for  Van  der  Pol 
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TABLE  n 

Certified  normalized  volume  on  examples  SI -S6  from  [4].  The  vector  fields  for  Examples  S3  and  S4  have 

DEGREE  EQUAL  TO  5,  WHILE  ALL  OTHERS  HAVE  DEGREE  EQUAL  TO  3. 
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2 

- 

TABLE  m 


Computation  statistics  for  the  benchmark  example 
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successes 

worst  (in  300) 

over  100 

time  (sec) 

2 

13 
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0.70 
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25 
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1.12 

4 

48 
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0.99999 

2.14 

6 
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297 

0.99997 

0.99998 

11.2 

8 

420 

300 

0.99989 

0.99992 

99.7 

Figure  Captions 


Fig.  1.  Provable  ROA  using  pointwise  maximum  of  two  polynomial  functions,  with  shape  factor 
xTR\X 

Fig.  2.  Provable  ROA,  from  [11],  [7]  and  [8]. 

Fig.  3.  VDP;  Level  sets  of  two  6th  degree  polynomials  at  Vi,  =  1 

Fig.  4.  Provable  ROA  using  pointwise  max  of  two  polynomial  functions,  with  shape  factor  x1  R^x 

Fig.  5.  Trajectories  and  level  sets  V  <  1.  LEFT  panel:  example  from  Section  III-C;  RIGHT 
panel:  example  from  Section  III-D 
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Local  Robust  Performance  Analysis  for  Nonlinear  Dynamical  Systems 


Ufuk  Topcu  and  Andrew  Packard 


Abstract — We  propose  a  computational  method  for  local 
robust  performance  analysis  of  nonlinear  systems  with  poly¬ 
nomial  dynamics.  Specifically,  we  characterize  upper  bounds 
for  local  £2  — ♦  £2  input-output  gains  using  polynomial  Lya¬ 
punov/storage  functions  satisfying  certain  dissipation  inequali¬ 
ties  and  compute  safe  approximations  for  these  upper  bounds 
via  sum-of-squares  programming  problems.  We  consider  both 
bounded  parametric  uncertainties  and  bounded  uncertainties 
due  to  unmodeled  dynamics. 

I.  Introduction 

We  consider  the  problem  of  quantifying  robust  perfor¬ 
mance  properties  of  uncertain  nonlinear  dynamical  systems 
with  polynomial  vector  fields  around  asymptotically  stable 
equilibrium  points.  The  amount  of  amplification  of  bounded 
£2  input  norms  at  the  output  channels  is  used  as  a  measure 
of  performance.  Two  types  of  uncertainties  are  considered: 

(1)  bounded  uncertainties  due  to  unmodeled  dynamics  and 

(2)  bounded  parametric  uncertainties.  Following  [1],  [2],  [3], 
we  characterize  upper  bounds  on  local  input-output  gains 
due  to  bounded  £2  disturbances  by  Lyapunov/storage  func¬ 
tions  which  satisfy  certain  “local”  dissipation  inequalities 
[4].  Similar  problems  were  studied  in  [1],  [2],  [5],  [6], 

[7]  mainly  for  systems  with  no  uncertainty.  Input-output 
properties  of  uncertain  nonlinear  systems  were  examined  in 

[8]  (for  discrete  time  nonlinear  systems  with  a  finite-time 
horizon  performance  metric)  and  [9]  (input-output  gains  for 
sufficiently  small  input  signals). 

In  this  paper,  we  use  polynomial  Lyapunov/storage  func¬ 
tion  candidates,  simple  generalizations  of  the  S-procedure 
[10],  and  sum-of-squares  (SOS)  relaxations  for  polynomial 
nonnegativity  [11]  and  compute  upper  bounds  on  the  input- 
output  gains  via  (bilinear)  SOS  programming  problems. 
Uncertainties  due  to  unmodeled  dynamics  are  accounted  for 
in  the  setting  [12]  shown  in  Figure  1  where  M  models 
the  nominal  part  and  $  is  an  unknown  operator  satisfying 
certain  relations  between  the  input  2  and  the  output  W2-  The 
objective  is  to  compute  upper  bounds  on  the  £2  norm  of  the 
exogenous  output  e  in  terms  of  the  £2  norm  of  the  exogenous 
input  w\.  The  approach  is  composed  of  two  steps:  first  bound 
the  £2  norm  of  the  internal  input  W2  to  M  in  terms  of  the  £2 
norm  of  w\  and  then  perform  an  input-output  gain  analysis 
for  M  from  the  inputs  (wi,W2)  to  the  output  e. 

The  approach  for  the  bounded  parametric  uncertainties  is 
similar  to  that  developed  in  [13],  [14]  in  the  context  of 

U  Topcu  is  with  Control  and  Dynamical  Systems  at  California  Institute 
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Fig.  1.  Input-output  system  with  the  feedback  interconnection  of  <£  and 
Af. 


robust  region-of-attraction  analysis.  Namely,  a  parameter- 
independent  Lyapunov/storage  function  is  used  to  charac¬ 
terize  input-output  properties  of  uncertain  systems  over  the 
entire  set  of  admissible  values  of  uncertain  parameters.  The 
input-output  relations  characterized  by  a  single  parameter- 
independent  certificate  may  be  more  conservative  compared 
to  those  by  parameter-dependent  certificates.  This  potential 
conservatism  is  simply  reduced  by  partitioning  the  set  of  un¬ 
certain  parameters  into  subregions  and  computing  parameter- 
independent  certificates  for  each  subregion.  The  partition  of 
the  uncertainty  set  can  be  refined  following  ideas  parallel  to 
branch-and-bound  algorithm  [15]  to  further  reduce  the  con¬ 
servatism.  Although  it  is  simplistic  (compared  to  techniques 
based  on  parameter-dependent  Lyapunov  functions),  this  ap¬ 
proach  offers  certain  computational  advantages  as  discussed 
in  [14]  for  robust  region-of-attraction  analysis.  In  fact,  in 
robustness  analysis  involving  time-invariant  unknown  param¬ 
eters,  it  is  common,  [16],  [17],  to  combine  easily-computable 
sufficient  conditions  with  branch-and-bound  strategies,  often 
yielding  improved  analysis  results. 

The  rest  of  the  paper  is  organized  as  follows:  A  character¬ 
ization  of  upper  bounds  for  £2  — ►  £2  input-output  gains 
by  Lyapunov/storage  functions  is  discussed  in  section  II. 
Section  III  is  devoted  to  the  development  of  the  results  for  the 
case  with  unmodeled  dynamics  and  this  is  followed  by  the 
method  to  account  for  parametric  uncertainties  in  section  IV. 
Implementation  details  are  given  in  section  V.  Demonstration 
of  the  methodology  with  examples  in  section  VI  precedes  the 
concluding  remarks. 

Notation:  For  £  £  7Zn ,  £  y  0  means  that  &  >  0  for 
&  =  !,*•  ,  n.  For  Q  =  QT  €  TZnxny  Q  0  (Q  >-  0)  means 
that  £rQ£  >  0  (>  0)  for  all  £  c  7 ZN .  R[£]  represents  the  set 
of  polynomials  in  £  with  real  coefficients.  The  subset  £[£]  := 

{7T  =  7[f  +  7 rf  + - 1-  :  7T1 ,  •  ■  •  ,  7TM  G  R[£]}  of  R[£] 

is  the  set  of  SOS  polynomials.  For  ~  £  R[f],  d{~)  denotes 
the  degree  of  n.  For  77  >  0  and  a  function  g  :  7Zn  — >  7£, 
define  the  77-sublevel  set  of  g  as 

:=  {x  £  Rn  :  g(x)  <  rj}. 


In  several  places,  a  relationship  between  an  algebraic  con¬ 
dition  on  some  real  variables  and  state  properties  of  a  dy¬ 
namical  system  is  claimed,  and  same  symbol  for  a  particular 
real  variable  in  the  algebraic  statement  as  well  as  the  state 
of  the  dynamical  system  is  used.  This  could  be  a  source  of 
confusion,  so  care  on  the  reader’s  part  is  required.  < 


II.  Upper  bounds  on  the  £2  — ►  £2  input-output 

GAIN 


Consider  the  dynamical  system  governed  by 

x(t)  =  f(x(t),w(t)) 
z(t)  =  h(x(t)), 


where  x(£)  G  7£n,  G  TV1™ ,  and  /  is  a  n-vector  with 
elements  in  R[(x,u;)]  such  that  /(0,0)  =  0  and  h  is  an 
n2-vector  with  elements  in  R[x]  such  that  /i(0)  =  0.  Let 
<f>(t\ xo,  w)  denote  the  solution  to  (1 )  at  time  t  with  the  initial 
condition  x(0)  =  xo  driven  by  the  input/disturbance  w.  For 
a  piecewise  continuous  map  u  :  [0,  00)  — ►  7£m,  define  the 
(truncated)  £2  norm  as 


IMkr  := 


u(t)Tu(t)dt. 


For  notational  simplicity,  denote  ||u||2i00  by  ||u||2. 


A.  Upper  bounds  on  local  £2  — ►  £2  gain 

Lemma  II.l.  [2]  If  there  exist  a  real  scalar  7  >  0  and  a 
continuously  differentiable  function  V  such  that ,  for  R  >  0, 

2  is  bounded , 

V'(O)  =  0  and  V (x)  >  0  for  all  nonzero  x  G  72n,  (2) 

W/(x,  w)  <  wTw  —  7 ~2zT  z  Vx  G  fVi?2  an/d  Vu;  G  1Znw , 

(3) 

then  it  holds  that  for  the  system  in  (1)  and  for  all  T  >  0 
IM|2,t  <  R  and  x(0)  =  0  =►  \\z\\2,t  <  TlMkr-  (4) 


< 


In  other  words,  7  is  a  local  upper  bound  for  the  input- 
output  gain  for  the  system  in  (1).  We  call  7  to  be  a  local 
upper  bound  because  the  upper  bound  on  the  norm  of  the 
output  z  is  only  supposed  to  hold  whenever  the  norm  of  the 
input  is  bounded  by  R.  This  is  unlike  the  input-output  gains 
for  linear  systems  which  hold  for  all  values  of  input  norms. 

Let  7  >  0  be  fixed  and  V  be  the  space  of  continuously 
differentiable  functions.  Define  Rc7,opt(Vi'lf)  be  the  max¬ 
imum  value  of  R  such  that  the  conditions  in  Lemma  II.l 
hold  for  some  V  G  V.  Let  Vpo/y  be  a  subset  of  V  that 
is  composed  of  all  polynomials  in  x  of  some  fixed  finite 
degree  (omitted  in  notation).  By  restricting  the  search  for  V 
satisfying  the  conditions  in  Lemma  II.l  to  Vpoiy,  utilizing 
a  generalization  of  the  S-procedure  (see  Lemma  VUI.l  in 
the  Appendix)  to  obtain  sufficient  conditions  for  the  set  con¬ 
tainment  constraints  in  Lemma  II.l  and  SOS  relaxations  for 
polynomial  nonnegativity,  the  following  proposition  provides 
an  upper  bound  on  Rc7,opt(V,  7)- 


Proposition  II.L  For  given  7  >  0  and  positive  definite 
polynomial  l  let  Rc2  be  defined  through 


R2cSVpoiy,S,l,l)  :=  max 

2  V£Vpolv,R>0,s£S 


R2 


s.  t. 


V (0)  =  0,  s  G  E[(x,tu)], 

v-ie  E[x], 

-  [(R2  -  V)s  4-  W /  -  wTw  +  7 ~2zTz\ 


(5) 

(6) 
(7) 


G  E[(x,  w)]. 


where  Vpoiy  C  V  is  as  defined  above  and  S  is  a 
prescribed  finite-dimensional  subset  of  R[(x,w)].  Them 
Rc2(Vpoly,S,7il)  <  #£2i0pt(V,7).  < 

Note  that  Rc7,opt  depends  on  7  and  V  and  Rc2  depends 
on  Vpoiy 1  <$,  e,  U  and  7.  Hereafter,  this  dependence  will  not 
be  notated  explicitly  unless  it  causes  confusion. 

The  optimization  problem  in  Proposition  II.l  can  be  cast  in 
a  bilinear  SDP  (i.e.,  nonconvex  in  general).  Bilinear  SDPs  are 
known  to  be  harder  than  linear  ones  [18].  Consequently,  the 
state-of-the-art  of  the  solvers  for  bilinear  SDPs  is  far  behind 
that  for  the  linear  ones  and  methods  for  bilinear  SDPs  are 
generally  based  on  heuristics  such  as  coordinate- wise  affine 
search  or  specialized  solvers  e.g.  PENBMI[19].  Although 
these  techniques  are  local  search  schemes  and  convergence  to 
a  global  optimum  is  not  guaranteed,  coupled  with  efficient 
initializations,  they  have  been  effectively  used  for  several 
system  analysis  questions  [20],  [21].  For  the  examples  in 
this  paper,  we  use  a  coordinatewise  affine  search  scheme  as 
detailed  in  section  V. 

For  given  7  >  0,  the  optimization  problem  in  Propo¬ 
sition  II.l  maximizes  R  (that  can  be  verified  through  the 
families  of  admissible  Lyapunov  function  candidates  (U) 
and  S-procedure  multipliers  (s))  such  that  ||z||2  <  7||u;||2 
whenever  ||w||2  <  R .  One  can  also  choose  to  minimize  7 
for  a  given  value  of  R  and  this  can  be  formulated  as  an 
optimization  problem  similar  to  that  in  Proposition  II.l  with 
minor  changes. 

III.  Robust  performance  in  the  presence  of 

UNMODELED  DYNAMICS 

Consider  the  input-output  system  in  Figure  1.  Let 

x{t)  =  f(x(t),Wi(t),W2(t)) 

z(t)  =hl(x(t))  (14) 

e(t)  =  h2{x(t)) 

be  a  realization  of  M,  where  W\  ( t )  G  ,  w2(f)  G  , 
/  is  an  n-vector  of  polynomials  in  R[(x, tui,u;2)]  with 
/( 0, 0, 0)  =  0,  h\  and  h2  are  nz  and  ne  dimensional  vectors 
with  entries  in  R[x]  satisfying  h\(0)  =  0  and  h2(0)  =  0. 
Furthermore,  assume  that  is  causal  and,  starting  from  rest, 
satisfies 

||3Hz)||2,T  =  ||u>2||2,T  <  IM|2,T  (15) 


for  all  T  >  0.  Lemma  III.l  provides  a  bound  on  ||tu2||2,T  in 
terms  of  ||u?i||2,t  for  T  >  0.  In  the  following  proposition, 
this  result  will  be  used  to  establish  a  local  upper  bound  on 
the  norm  of  the  exogenous  output  e  in  terms  of  the  norm  of 
the  exogenous  input  w\ . 


R?  :=  max 

vevpoiv,R>o,se5 


R2  subject  to 


V(0)  =  0,  s  G  E((x,t0i,tU2)], 

-  [(R2  -  V)s  +  VVf(x,w1,wz)  ~  (02w\wi  +wZw2)  +  a-2zTz]  G  E[(x,ti>i,ti;2)], 


(8) 

(9) 

(10) 


7i  :=  min 

QeVpoiv,7>0,s65 


7  subject  to 
<5(0)  =0,  seE[(i,u)i,w2)],  <5-(gE(x], 

-  [(l^a7  -Q)  s  +  'VQf(x,w  1,^2)  -  (02w[wi  +W2W2)  +7-2ere 


€  £[(x,u’i,u;2)]. 


(ID 

(12) 

(13) 


Lemma  IIL1.  For  R  >  0,  0  <  c*  <  1  and  f3  >  0,  if  there 
exists  a  continuously  differentiable ,  positive  definite  function 
V  such  that  V^O)  =  0,  ^lv,R2  w  bounded ,  and 

Wf(x,  w\,  w2)  <  02wjw\  +  w^w2  — \zT z 

Or 

for  all  x  e  w\  £  ,  and  w2  G  7lnuii ,  then  for  <£ 

starting  from  rest  and  for  all  T  >0 

x(0)  =  0  and  ||u;i||2,r  <  =*  11^2  II 2, r  <  ,  *  • 

P  V 1  —  or 

< 


Proof:  While  solutions  to  (14)  exist,  for  T  >  0 

02\\wl\\lyT  +  ll^lla.T  “  J^lklla,T 

<  p2\\wi ||2tr  “  II 2 II 2, r  —  P2\\wi lli.r- 


(i6) 


Since  f V,i?2  is  bounded,  as  long  as 


ll^i  lb,T  < 


v  r2  -  v{xm 


solutions  to  (14)  exist  for  all  T  >  0  and  satisfy 


a2/?2 


a 


IMIl,r  <  M&r  <  Mir  + 

In  particular,  for  x(0)  =  0,  V(x(0))  =  and  ||u’2||«;,t  < 

□ 

Proposition  III.l.  In  addition  to  the  conditions  in  Lemma 
///./,  if  there  exists  a  continuously  differentiable,  positive 
definite  function  Q  such  that  Q( 0)  =  0  and 

VQf{x,wi,w2)  <  (32wjwi  +  W2  u'2  —  ^ cT6 

for  all  x  G  /(i_a2)>  u>i  €  ,  w2  G  7 V1™* ,  then  for 

&  starting  from  rest  and  for  all  T  >  0 

||tUi||2,T  <  ^  and  x(0)  =  0  =>  ||e|| 2,r  <  jR/Vl  —  a2. 


Proof:  By  Lemma  III.l, 

r2 

||^i||2,r  <  H//3  =>  /32||^i||2,t  +  Ilu’2|l2,r  ^  \  _  'a2' 
Consequently,  the  result  follows  from  Lemma  D.  1.  □ 


Lemma  III.l  and  Proposition  III.l  can  be  used  to  construct 
relations  between  ||u>i||2,r  and  ||e||2,r-  Similar  to  Propo¬ 
sition  II.  1,  one  can  obtain  sufficient  conditions  for  those 
in  Lemma  III.l  and  Proposition  III.  1  using  Lemma  VIII.  1 
and  SOS  relaxations  for  polynomial  nonnegativity.  For  given 
/?  >  0,  a  >  0,  and  l(x)  =  txrx  (  with  e  >  0  fixed),  solve 
the  problems  in  (8)-(10)  and  (11)-(13).  Then,  for  <£  starting 
from  rest  and  for  all  T  >  0 

x(0)  =  0  and  Ikxb.r  <  ^  =*  ||e||2,T  <  -7=*. 

P  V  1  —  Cl1 


When  is  unknown  but  a  global  gain  relation  between  its 
inputs  and  outputs  is  known,  then  the  results  of  this  section 
provide  a  framework  for  robust  performance  analysis  for  the 
feedback  interconnection  between  M  and  <£.  On  the  other 
hand,  even  when  the  operator  is  known,  the  procedure 
outlined  in  this  section  can  be  used  as  a  framework  for 
compositional  performance  analysis.  Note  that  conditions 
in  Lemma  III.l  and  Proposition  III.l  do  not  involve  the 
states  of  (the  realization  of)  <£.  When  $  has  the  state  space 
realization  x2(t)  =  f2{x2{t),  z{t)),  then  the  gain  relation 
||$(z)||2,r  <  IMl2,r  can  be  established  by  determining 
a  positive  definite,  continuously  differentiable  function  V2 
satisfying 

w2/2(x2,  z)  <  ZTZ  -  wlw2  Vx2  G  nn\  Vz  G  TZn- 

which  does  not  depend  on  the  states  x\  of  M.  Consequently, 
Lemma  III.l  and  Proposition  III.l  enable  performance  anal¬ 
ysis  for  the  feedback  interconnection  between  M  and  $ 
based  on  the  input-output  properties  of  individual  blocks. 
Of  course,  this  analysis  may  be  conservative.  Nevertheless, 
compositional  analysis  may  be  a  fruitful  direction  which 
extends  the  applicability  of  SOS  programming  based  non¬ 
linear  analysis  tools  for  reasonably  larger  dimensional  sys¬ 
tems  whenever  it  is  possible  to  establish  an  interconnection 
structure  as  in  Figure  1.  Furthermore,  it  may  be  possible 
to  refine  the  input-output  relations  between  W\  and  e  by 
using  transformations  at  the  interconnections  similar  to  the 
D  scales  in  linear  robustness  analysis  [12]. 


IV.  Robust  performance  analysis  in  the  presence 

OF  PARAMETRIC  UNCERTAINTIES 

We  now  generalize  the  development  in  section  II  to  the 
case  where  the  vector  field  contains  unknown  but  fixed  and 
bounded  parameters.  Following  the  methodology  proposed 
in  [14]  in  the  context  of  robust  stability  analysis,  we  first 
restrict  our  attention  to 

x(t)  =f(x,w,S) 

:=  f0(x(t),  w(t))  +  ,  6ifi(x(t),  w{t))  (17) 

z(t)  =  h(x(t)), 

where  /o ,  /i ,  •  • . ,  /m  are  n- vectors  with  elements  in 
R[(x,tu)]  such  that  /o(0,0, 8)  =  /i(0,0,5)  =  ...  = 
/m(0,0, 8)  =  0,  for  all  8  £  A  C  Tlm,  and  A  is  a  known 
bounded  polytope.  Let  0(£;xo,tu,5)  denote  the  solution  of 
(17)  for  8  at  time  t  with  the  initial  condition  x(0)  =  xo 
driven  by  the  input/disturbance  w  and  £a  denote  the  set  of 
vertices  of  A. 

Proposition  IV.l.  If  there  exist  a  real  scalar  7  >  0  and  a 
continuously  differentiable  function  V  such  that  V  (0)  =  0, 
V(x)  >  0  for  all  nonzero  x  £  1Zn,  is  bounded,  and 

W  f(x,  w,  5)  <  wTw  —  ^~2zT  z  (18) 

for  all  x  £  w  £  Rnw,  and  8  £  £&,  then  the  system 

in  (17)  with  x(0)  =  0  satisfies  \\z\\2  <  7||w||2  whenever 
||tu||2  <  R  and  8  £  A,  < 

Proof:  Since  the  vector  field  is  affine  in  8  and 

A  is  a  bounded  polytope,  it  follows  that,  for  <5  £  A, 
W f(x,w,8)  <  wTw  for  all  x  £  Qy^R2,  w  €  ■  By 

Lemma  II. 1,  for  each  8  £  A,  ||z||2  =  |j/i(0(s 0, tx;, 5))||2  < 
7IHI2  whenever  ||tu||2  <  R-  □ 

Note  that  restricting  the  attention  to  affine  uncertainty  8 
dependence  and  polytopic  A,  Proposition  IV.l  enables  to 
compute  upper  bounds  on  £2  —♦  £2  gain  for  the  system 
in  (17)  by  imposing  the  conditions  in  (18)  at  finitely  many 
8  £  £  a  instead  of  at  infinitely  many  8  £  A.  Furthermore, 
sufficient  conditions  for  those  in  Proposition  IV.l  can  be 
obtained  using  Lemma  Vffl.I  and  SOS  relaxations. 

The  approach  proposed  here  is  restrictive:  (1)  only  affine 
dependence  on  8  and  polytopic  A  are  allowed  (2)  SOS 
relaxations  for  the  conditions  in  Proposition  IV.l  may  include 
a  large  number  of  SDP  constraints  (3)  single  (^-independent) 
Lyapunov/storage  function  is  to  certify  properties  for  an 
entire  family  of  systems.  These  limitations  can  be  partially 
alleviated  using  techniques  proposed  in  [14]  in  the  context  of 
robust  region-of-attraction  analysis.  For  example,  polynomial 
dependence  on  8  in  the  vector  field  and  the  output  map  can 
be  handled  by  covering  the  graph  of  non-affine  functions  8 
(in  the  conditions  in  Proposition  IV.l)  by  bounded  polytopes. 
Furthermore,  the  fact  that  constraints  in  the  SOS  relaxations 
for  the  conditions  in  IV.l  are  only  coupled  through  the 
Lyapunov/storage  functions  (which  include  relatively  small 
portion  of  all  decision  variables  in  associated  SDPs)  can 
be  exploited  through  a  suboptimal  two-step  procedure:  pick 
a  point  in  A,  compute  a  Lyapunov/storage  function  for 
the  system  corresponding  to  that  point,  and  then  in  the 


second  step  determine  an  input-output  relation  certified  by 
the  Lyapunov/storage  function  determined  (fixed)  in  the 
first  step  which  holds  for  the  entire  family  of  admissi¬ 
ble  systems.  This  procedure  effectively  decouples  the  large 
number  of  constraints  in  the  second  step  enabling  use  of 
trivial  parallelization.  Finally,  conservatism  (due  to  using  a 
single  parameter-independent  Lyapunov  function  and  due  to 
the  suboptimal  two-step  procedure)  can  be  reduced  by  an 
informal  branch- and-bound  type  refinement  procedure  where 
A  is  partitioned  into  smaller  subregions  and  a  different 
Lyapunov/storage  function  is  computed  for  each  subregion. 
See  [14]  which  develops  a  similar  methodology  in  the  context 
of  robust  region-of-attraction  analysis. 

V.  Implementation  Issues 

The  SOS  relaxations  in  (5)-(7)  lead  to  bilinear  SDPs  due 
to  the  multiplication  between  the  decision  variables  in  V 
and  the  multipliers.  Therefore,  solution  techniques  for  these 
problems  are  usually  limited  to  local  search  schemes  such 
as  PENBMI  [19]  or  coordinate- wise  affine  search  based 
on  the  observation  that,  for  given  V  and  i?,  constraints 
in  these  problems  are  affine  in  the  decision  variables  in 
the  multipliers.  For  example,  one  can  obtain  a  suboptimal 
solution  for  the  problem  in  (5)-(7)  by  altematingly  solving 
the  following  two  problems  until  a  maximum  number  of 
iterations  is  reacaed  or  the  increase  in  the  value  of  certified 
R  becomes  smaller  than  a  pre-specified  tolerance  .  For  given 
V 

max  R2  subject  to  s  £  T,{(x,w)}  and  (7),  (19) 

K>0,s€S 

which  can  be  solved  using  an  off-the-shelf  affine  SDP  solver 
through  a  line  search  on  /?,  and  for  given  (feasible)  multiplier 
s 

max  R2  subject  to  V -l  £  E[x),  V(0)  =  0,  and  (7). 

R>0y  €Vp0jv 

(20) 

Furthermore,  by  a  change  of  variables,  it  is  possible  to 
iterate  without  a  line  search  in  the  first  step.  Indeed,  for 
/?  >  0,  if  the  problem  in  (5)-(7)  has  the  solution  Vj  and 
sj,  then 

i  =  min  -^7  subject  to 

/C€Vpoivil/*2>0,S€S  Rl 

s  e  £[(x,u;)],  K  e  R[x],  K  -  h/F?  e  £[x], 

-  [(1  —  K)s  +  VKf  -  -gj(wTw  -  7 -2zTz)] 

e  £[(x,u’)]. 

(21) 

Note  that  for  given  K  constraints  in  (21)  are  affine  in  1  /R  ' 
and  s.  In  fact,  optimal  values  of  s  and  K  are  s  =  Rjsi,  and 
K  =  Vi /Rl 

VI  Examples 

Consider  the  controlled  short  period  aircraft  dynamics  in 
Figure  2  where  xp  :=  [xi  X2  X3]r,  Xi,  X2,  and  X3  denote 
the  pitch  rate,  the  angle  of  attack,  and  the  pitch  angle. 


Fig.  2.  Controlled  short  period  aircraft  dynamics  with  unmodeled  dynamics. 


respectively,  and 


c(xp) 

'  £[xp  +  b\ 

Xp  = 

q(xp) 

+ 

62 

xi 

0 

where,  c  and  q  are  cubic  and  quadratic  polynomials,  respec- 
tively,  ib  S  R3,  f>i  and  62  are  real  scalars  (see  [22]  for 
the  values  of  the  missing  parameters).  The  plant  output  is 
[xi  x3f  •  The  input  u  to  the  plant  is 

u  =  1.25t>  -f  W\  4-  V02 

where  v ,  the  elevator  deflection,  is  the  controller  output 
determined  by 

x4  =  — 0.864t/i  -  0.321y2 
v  =  2x4, 

where  X4  is  the  controller  state.  Assume  that  <£  :  1Z  — >  H 
satisfies,  starting  from  rest, 

||$(z)l|2, T  =  |k2||2,T  <  |kl|2,T 

for  all  T  >  0.  We  performed  the  following  analysis: 

(i)  For  several  values  of  a  €  [0.55,0.9],  solve  the  prob¬ 
lems  in  (8M10)  and  (11M13). 

(ii)  Apply  linearized  robust  performance  analysis  for  the 
feedback  interconnection  [23]  and  fit  a  first  order  stable 
minimum  phase  transfer  function,  say  H(s)y  to  the 
optimal  D-scales.  For  several  values  of  a  e  [0.55,0.9], 
solve  the  problems  in  (8)-(  10)  and  (11 )-( 1 3)  for  the 
system  HMH~l  with  a  minimal  realization  for  H. 

(iii)  Solve  the  problem  in  Proposition  II.  1  for  the  system 
with  no  uncertainty  for  several  values  of  7. 

Figure  3  shows  the  £2  norms  of  the  exogenous  outputs  e 
versus  the  £2  norms  of  the  exogenous  inputs  w\  in  each  of 
these  cases:  («)  with  marker  “+”,  («)  with  marker  and 
(m)  with  marker  “x”. 

Figure  3  illustrates  the  trade  off  between  the  robustness 
and  performance:  As  a  gets  larger,  the  gap  between  the 
nominal  performance  level  and  the  “robust”  performance 
level  increases  deduced  from  the  divergence  between  the 
curve  with  “+”  and  other  two  curves. 

VII.  Conclusion 

We  proposed  a  computational  method  for  local  robust 
performance  analysis  of  nonlinear  systems  with  polynomial 
dynamics.  Specifically,  we  characterized  upper  bounds 
for  local  £2  — >  £2  input-output  gains  using  polynomial 
Lyapunov/storage  functions  satisfying  certain  dissipation 


inequalities  and  computed  safe  approximations  for  these 
upper  bounds  via  sum-of-squares  programming  problems. 
We  considered  both  bounded  parametric  uncertainties  and 
bounded  uncertainties  due  to  unmodeled  dynamics. 
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VIII.  Appendix 

The  following  lemma  is  a  simple  generalization  of  the  S- 
procedure  [10]  and  is  used  to  obtain  sufficient  conditions  for 
certain  set  containment  constraints  throughout  the  paper. 

Lemma  VUI.l.  For  goy9iy---  >Pm  €  R[x],  if  there  exist 
s\y  •  •  •  ,  sm  €  E[x]  such  that 

m 

go  ~^2si9i  s  S[x], 

1=1 

then 

{x£Un  :gi(x),...,gm{x)  >  0} 

C{i£Rn:  g0{x)  >  0}  . 

<3 
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Linearized  Analysis  versus  Optimization-based  Nonlinear  Analysis  for 

Nonlinear  Systems 

Ufuk  Topcu  and  Andrew  Packard 


Abstract — For  autonomous  nonlinear  systems  stability  and 
input-output  properties  in  small  enough  (infinitesimally  small) 
neighborhoods  of  (linearly)  asymptotically  stable  equilibrium 
points  can  be  inferred  from  the  properties  of  the  linearized  dy¬ 
namics.  On  the  other  hand,  generalizations  of  the  S-procedure 
and  sum-of-squares  programming  promise  a  framework  poten¬ 
tially  capable  of  generating  certificates  valid  over  quantifiable, 
finite  size  neighborhoods  of  the  equilibrium  points.  However, 
this  procedure  involves  multiple  relaxations  (unidirectional 
implications).  Therefore,  it  is  not  obvious  if  the  sum-of- 
squares  programming  based  nonlinear  analysis  can  return  a 
feasible  answer  whenever  linearization  based  analysis  does. 
Here,  we  prove  that,  for  a  restricted  but  practically  useful 
class  of  systems,  conditions  in  sum-of-squares  programming 
based  region-of-attraction,  reachability,  and  input-output  gain 
analyses  are  feasible  whenever  linearization  based  analysis  is 
conclusive.  Besides  the  theoretical  interest,  such  results  may 
lead  to  computationally  less  demanding,  potentially  more  con¬ 
servative  nonlinear  (compared  to  direct  use  of  sum-of-squares 
formulations)  analysis  tools. 

I.  Introduction 

Internal  stability,  input- to-state,  and  input- to-output  prop¬ 
erties  of  dynamical  systems  are  commonly  analyzed  by 
constructing  Lyapunov/storage  functions  satisfying  certain 
conditions  (such  as  dissipation  inequalities)  [1],  [2],  [3], 
[4].  Generalizations  of  the  S-procedure  [5],  [4]  and  sum-of- 
squares  (SOS)  relaxations  for  polynomial  nonnegativity  [6] 
provide  a  framework  for  the  search  of  such  Lyapunov/storage 
functions  for  systems  with  polynomial  vector  fields  based  on 
(linear  or  bilinear)  semidefinite  programming  (SDP)  prob¬ 
lems  [7],  [8],  [9],  [10],  [11],  [12],  [13],  [14],  [16],  [17],  [18]. 

On  the  other  hand,  it  is  well  known  that  if  there  exist  Lya¬ 
punov/storage  functions  for  the  linearized  dynamics  (around 
an  asymptotically  stable  equilibrium  point)  then,  by  certain 
continuity  assumptions,  these  functions  (always)  serve  as 
Lyapunov/storage  functions  for  the  nonlinear  system  possibly 
only  locally,  i.e„  corresponding  Lyapunov  or  dissipation 
inequalities  only  hold  in  a  “sufficiently  small”  neighborhood 
of  the  equilibrium  point.  The  promise  of  SOS  programming 
based  nonlinear  analysis  is  that  it  may  be  possible  to  con¬ 
struct  Lyapunov/storage  functions  that  satisfy  the  Lyapunov 
or  dissipation  inequalities  not  only  in  a  “sufficiently  small” 
neighborhood  of  the  equilibrium  point  but  also  over  quan¬ 
tifiable,  non-trivial  subsets  of  the  state  space.  However,  the 
transformation  from  system  analysis  questions  to  correspond¬ 
ing  SDP  problems  (in  nonlinear  analysis)  involves  a  series 
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of  sufficient  (but  not  necessarily  necessary)  conditions.  For 
example,  except  certain  special  or  hypothetical  cases,  S- 
procedure  is  not  lossless  and  not  all  nonnegative  polynomials 
are  SOS  [9],  [6],  [19].  Therefore,  it  is  not  obvious  if  (SOS 
programming  based)  nonlinear  analysis  yields  a  certificate 
for  the  nonlinear  system  whenever  the  linear  analysis  does. 

In  this  paper,  we  propose  conditions  for  the  feasibility 
of  SDP  problems  (equivalently  SOS  programming  prob¬ 
lems),  proposed  in  [5],  [20],  [7]  in  the  context  of  stability 
robustness,  reachability,  and  input-output  gain  analyses  of 
nonlinear  systems  around  asymptotically  stable  equilibrium 
points,  based  on  the  properties  of  the  corresponding  lin¬ 
earized  dynamics.  We  focus  on  systems  with  cubic  polyno¬ 
mial  vectors  fields  mainly  due  to  practical  reasons.  Although 
SOS  programming  based  analysis  can  theoretically  be  used 
for  systems  with  polynomial  vector  fields  of  any  finite 
degree,  there  are  practical  bounds  on  the  degree  imposed  by 
the  capabilities  of  current  SDP  solvers  and  computational 
resources  (see  [7],  [14]  for  a  more  detailed  discussion). 
Therefore,  nonlinear  analysis  with  cubic  vectors  fields  is 
a  pragmatic  extension  for  linearization  based  analysis  with 
tighter  approximations  for  the  actual  dynamics  and  richer 
families  of  Lyapunov/storage  functions. 

The  motivation  is  primarily  theoretical,  showing  that  the 
optimization-based  (S-procedure/SOS)  methods  for  nonlin¬ 
ear  analysis  (as  proposed  in  [5],  [20],  [7])  always  involve 
feasible  bilinear  SDP  problems  whenever  the  linearization 
is  asymptotically  stable.  Furthermore,  these  results  may  also 
have  some  limited  practical  value  in  actually  constructing 
(possibly  conservative)  quantitative  results  for  the  nonlinear 
system  as  outlined  in  section  VI. 

Notation:  For  £  6  7£n,  £  X  0  means  that  £*  >  0  for 
k  =  I,---  , n.  For  Q  =  QT  €  Knxny  Q  fc  0  (Q  y  0) 
means  that  £TQ£  >  0  (>  0)  for  all  £  e  1Zn.  For  xi  e  1Zni 
and  X2  €  'RJ12,  [xi;x2]  €  7£ni+na  denotes  the  concatenation 
of  Xj  and  X2-  R[£]  represents  the  set  of  polynomials  in  £ 
with  real  coefficients.  The  subset  E[£]  :=  {7r  =  7Tj  +  -f 

- b  :  TTi,  •  •  *  ,  km  6  R[£]}  of  R[£]  is  the  set  of  SOS 

polynomials.  For  rj  >  0  and  a  function  g  :  RJ1  — ►  7 Z,  define 
the  rj-sublevel  set  of  g  as 

:=  {x  £Hn  :  g(x)  <  rj}. 

For  a  piecewise  continuous  map  u  :  [0,  oo)  — ►  7 Zm,  define 
the  £2  norm  as 

Ml.  :=  J l  u{t)Tu{t)dt. 


In  several  places,  a  relationship  between  an  algebraic  con¬ 
dition  on  some  real  variables  and  state  properties  of  a  dy¬ 
namical  system  is  claimed,  and  same  symbol  for  a  particular 
real  variable  in  the  algebraic  statement  as  well  as  the  state 
of  the  dynamical  system  is  used.  This  could  be  a  source  of 
confusion,  so  care  on  the  reader’s  part  is  required.  < 


for  t  =  l , ,n  and  Xiyj  for  i  =  1, . .  ,n  and  j  =  1, . . .  m 
with  no  repetition,  <\ 

Although  z  (as  defined  above)depends  on  x  and/or  y, 
this  dependence  will  not  be  explicitly  notated  whenever  it 
is  convenient  and  does  not  cause  confusion. 


II.  Preliminaries 


III.  C2  ->  C2  INPUT-OUTPUT  GAIN  ANALYSIS 


Following  two  lemmas  are  straightforward  generalizations 
of  the  S-procedure  [4].  See  [21],  [7]  for  the  proofs. 

Lemma  ILL  Given  yo>  <7i»  *  *  *  i  <7m  €  R[x],  if  there 

exist  Si,--  -  ,sm  €  £[x]  such  that  g0  —  YliL isi9*  € 
£[x],  then  {x£lZn  :  gi(x), . . .  ,ym(x)  >  0} 

C  {x  e7ln  :  gQ(x)  >  0}  . 

Lemma  II.2.  Given  yoi  9\,  92  €  R[x]  such  that  go  is  positive 
definite  and  yo(0)  =  0,  if  there  exist  $\,s2  €  £[x]  such  that 
9iSi  “b 92$2  -90  €  £[x],  then  {x  €  TZn  :  gi(x)  <  0}\{0} 
C  {x  €  TZn  :  g2{x)  >  0}.  < 


The  following  fact  will  be  used  in  the  subsequent  sections. 

Lemma  1 13.  Let  Q  =  QT  €  lZnxn  be  positive  definite , 
/  :  lZn  — ►  1Z  be  defined  as  f(x)  =  xTQx,  Ci,...,Cm 
be  positive  real  numbers ,  and  g  :  7Zm  — >  1Z  be  defined  as 
g(y)  =  ciy\  4-  c2yl  +  . . .  +  Cmy^.  Then,  f{x)g(y)  can  be 
written  in  the  form 

f(x)g(y)  =  z(x,y)T  Hz(x,y), 

where  z(x,  y)  =  y  ®  x  and  H  >~  0.  < 


Proof: 

f(x)g(y) 


=  xTQx{ciyl  +  •  •  •  +  Omym) 

=  Y.?=\<*(yix)TQ(yix) 

=  z(x,y)T  Hz(x,y), 


where  H  =  HT  €  ftnmxnm  is 

c\Q 

H  :=  ' 


CnQ 


Clearly,  H  is  positive  definite  since  Q  is  positive  definite.  □ 

Lemma  IL4.  Let  Q  and  f  be  as  in  Lemma  113,  C\ , . . . ,  Cn 
be  positive  real  numbers,  and  g  :  7Zn  1Z  be  defined  as 
g(x)  =  Cix\  +  . . .  +  Cnx\.  Then,  f(x)g(x)  can  be  written 
in  the  form 

f(x)g(x)  =  z(x)T  Hz(x), 


where  z(x)  is  a  vector  of  all  monomials  of  the  form  Xiyj  for 
i  =  1, . . . ,  n  and  j  >  i  with  no  repetition,  < 

Lemma  II.5.  Let  Q  and  f  be  as  in  Lemma  IL3,  C\ , . . . ,  c„+m 
be  positive  real  numbers,  and  g  :  7£m+n  —*  1Z  be  defined  as 
g(x,y)  =  Ci2/i+c2j/2+-  •  •+cmy^J+cm+ixf+. .  .+cm+ni^. 
Then,  /(x)y(x,y)  can  be  written  in  the  form 

/(x)y(x,y)  =  z(x,  [x;  y])f  Hz(x,  [x;y]), 

where  z(x,  [x;  y})  is  a  vector  of  all  monomials  of  the  form  x2 


Consider  the  dynamical  system  governed  by 

x(t)  =  f(x(t),w(t)) 
y(t)  =  h(x(t )), 

where  x(t)  €  lZn ,  w(t)  €  7Znu\  and  /  is  a  n-vector  with 
elements  in  R[(x,tu)]  such  that  /( 0,0)  =  0  and  h  is  an 
ny-vector  with  elements  in  R[x]  such  that  /i(0)  =  0.  Let 
0(f;  xo,  w)  denote  the  solution  to  (1)  at  time  t  with  the  initial 
condition  x(0)  =  xo  driven  by  the  input/disturbance  w. 

Lemma  HL1.  [22]  If  there  exist  real  scalars  7  >  0  and 
R  >  0  and  a  continuously  differentiable  function  V  such 
that 

V (0)  =  0  and  V(x)  >  0  for  all  nonzero  x  €  lZn ,  (2) 

2  is  bounded ,  (3) 

Wf(x,w)  <  wTw  — 7“2yTy  Vxcfi^R2?  w£lZnrv, 

(4) 

then  for  the  system  in  (1) 

x(0)  =  0  and  ||tu||2  <  R  ||y||2  <  7|M|2-  (5) 

< 

In  other  words,  7  is  an  upper  bound  for  the  “local”  input- 
output  gain  for  the  system  in  (1).  For  given  7  >  0,  we 
restrict  V  to  be  a  polynomial  of  some  fixed  degree  and  use 
Proposition  m.l  to  compute  lower  bounds  on  the  maximum 
value  of  R  such  that  (2)-(4)  hold. 

Proposition  III.1.  [21]  For  given  7  >  0  and  positive  definite 
polynomial  1 1  satisfying  1 1(0)  =  0,  let  Rc7  be  defined 
through 

R2r  :=  max  R2  subject  to  (6) 

v'evpOiV,.R>0,s€«s 

1^(0)  =0,  8ieH[(x,w)l  (7) 

V-h  €E[x]f  (8) 

—  [(i?2  -  V)sx  +  W/(x,tc)  -wTw  +  7 ~2yTy]  (9) 

€  £[(x,u;)], 

where  Vpoiy  C  V  and  S  are  prescribed  finite -dimensional 
subsets  of  R[x].  Then, 

x(0)  =  0  and  \\w\\2  <  Rc2  ||y||2  <  711^112- 


Now,  consider  the  case  where  /  and  h  are  of  the  form 
f(x,w)  =  Ax +  Bw  +  f2(x)  -f  /3(x) 

+  (9i(x)  +  9z(x))Wi  (10) 

h(x)  =  Cx  -f-  /i2(x), 


< 


where  /2,  /3,  gu  <72,  and  h2  are  matrices  (of  appropriate 
dimension)  of  (purely)  quadratic,  cubic,  linear,  quadratic, 
and  quadratic  polynomials  in  their  arguments  respectively 
and  A,  B ,  and  C  are  matrices  (of  reals)  of  appropriate 
dimension.  Then,  the  following  proposition  gives  conditions 
on  the  feasibility  of  the  constraints  in  (6)-(9)  based  on  the 
analysis  of  the  corresponding  linearized  dynamics. 

Proposition  IIL2.  For  given  7  >  0,  l\(x)  =  xTL\X  with 
L\  >-  0,  /  and  h  in  the  form  in  (10),  if  there  exist  a  symmetric 
matrix  Q  and  e  >  0  such  that  Q  >z  L\  and 

[  AtQ  +  QA  +  1-*CtC  qb  1 
D°  -  [  BTQ  -I  J  -  ~a' 

then  the  constraints  in  (6)-(9)  are  feasible.  < 

Proof:  Define  V(x)  :=  xTQx.  Let  z  =  z(x,  [x;iu])  be 
as  defined  in  section  II.  Then,  there  exist  H  0,  Mu  and 
M2  such  that 

zT  Hz  =  ( wTw  4*  xTx)(xTQx) 
xTM\z  =  xTQ(/2(x)  +  gi(x)w)  +  xTCTh2(x) 
ztM2z  =  2xTQ(f$(x)  4-  g2(x)w)  4-  h2(x)Th2(x). 


Let  a  >  0  be  such  that 


Di  := 


—Do  —M\ 
-Mi  olH  -  M2 


hel 


and  R  y/ej(2ctj.  Define 

$i(x,w)  :=  a(xTx  4-  ivTw). 


Then,  V  —  l\  and  si  are  SOS.  Consider 

b(x,  w)  :=  —  [W/(x,  w)  —  wTw  4-  hT(x)h(x)] 
—  a(xTx  4*  wTw)  ( R 2  —  V)  , 

which  can  be  decomposed  as 

b(x,  w)  =  [x;  w\  z]tD2[x;  w;  z], 


where 


D2  :=  D 1  — 


'  aR2I 

0  ' 

h  el  - 

'  aR2I  0  ‘ 

0 

0 

0  0 

Hence,  b  is  SOS. 


*§'■ 


□ 


IV.  Reachability  analysis 

For  R  >  0  and  ||tt;||2  <  R,  the  set  Gr?  of  points  reachable 
from  the  origin  under  the  flow  of  (1)  is  defined  as 

gR2  :=  {4>(T',Q,w)  £lln  :  T>  0,  ||tu||2  <  R] . 

Lemma  IV.  1  adapted  from  a  Lyapunov-like  argument  in  [4, 
§6.1.1]  provides  a  characterization  of  sets  containing  Qri 
[5],  [22]. 

Lemma  IV.  1.  If  there  exists  a  continuously  differentiable 
function  V  such  that 

V(x)  >  0  for  all  x  €  7£n\{0}  with  V(0)  =  0,  (11) 

is  bounded,  (12) 

W/(x,u>)  <  wTw  Vx  G  nV'W,  w  G  nn'“,  (13) 


then  Gr?  Q  DVyR*  • 


For  given  /3  >  0  and  positive  definite,  convex  polynomial 
p,  the  following  proposition  provides  a  lower  bound  for  the 
maximum  value  of  R  such  that  Gr 2  C 


Proposition  IV.1.  [22]  Let  (3  >  0,  h  be  a  positive  definite 
polynomial  satisfying  lx( 0)  =  0,  Rreach  be  defined  through 


Rreach  ■■=  VPV  max  R  subject  to  (14) 

V€Vpoiv!^^0,$i€oi,32€o2 

V(0)  =0,  si  G  2[x],  and  S2  G  E[(x,io)],  (15) 

V-JiCEJ*],  (16) 

(/3-p)-(R*-V)st  GE[x],  (17) 

—  [(R2  —  V)s2  +  VV/(x,  w)  —  wTw]  G  E[(x, it’)], (18) 


where  Vpoiy  c  V  and  Si  are  prescribed  finite-dimensional 
subsets  of  R[x].  Then, 


GrI. 


Proposition  IV.2.  For  p(x)  =  xT Px  with  P  0,  /i(x)  = 
xT  L\X  with  L  i  >-  0,  am/  /  of  the  form  in  (10),  if  there  exist 
e  >  0  and  Q  >z  L\  such  that 


Do 


ATQ  +  QA  QB 
BtQ  -I 


then  the  constraints  in  (I4)-(18)  are  feasible.  < 


Proof:  Define  V(x)  :=  xTQx.  Let  z  =  z(x,  [x;  w})  be 
as  defined  in  section  II.  Then,  there  exist  H  >-  0  (by  Lemma 
II.3),  M\,  and  M2  such  that 

zT Hz  —  (wTw  4-  xTx)(xTQx) 
xT M\z  =  xTQ(f2(x)  4*  gi(x)w) 
zTM2z  =  2xTQ(/3(x)  4*  p2(x)ii?). 

Let  a  >  0  be  such  that 

D  •-  f  ~D°  1  >-  fT 

D'  —  [  -Ml  aH  -M2  J  -  7’ 

and  R  :=  y/e/(2a).  Define 

s\{x)  ;=  Amax  {D)/^min{Q) 

s2(x,w)  :=  a(xTx  4-  wTw). 

Then,  V  —  li,  si,  s2 ,  and  {/3  —p)  —  (R2  —  V*)si  are  SOS. 
Consider 

b(x,w)  :=  —VVf(x,w)+‘wTw—a(xTx+wTw)  ( R 2  -  V)  , 
which  can  be  decomposed  as 

b(x,w)  =  [x;  to;z]TD2[x;  w\  z], 


where 
D2  :=  D\ 


aR2I  0 

0  0 


he!  - 


aR2I  0 

0  0 


fciJ- 


Hence,  b  is  SOS. 


□ 


A.  Extensions  of  the  reachability  analysis  for  systems  with 
degenerate  linearization 


Consider  the  system 

x(t)  =  Amx(t)  4-  BAKj(t)x(t)  +  Ew(t) 
kx  =  —Txx(t)xT(t)PB, 

where  x(t)  £  7ZnJ  B  £  7Znxm,  w(t)  £  TZnX7luj ,  and  P,  E , 
A,  Am,  and  Tx  are  matrices  of  appropriate  dimension  with 
Hurwitz  Am.  The  dynamics  in  (19)  can  be  considered  as 
the  closed  loop  dynamics  for  the  system  x(£)  =  Amx{t)  + 
BAu(t)  regulated  to  the  origin  by  a  model  reference  adaptive 
controller  of  the  form[23] 

u(t)  =  kx(t)x(t) 

in  the  presence  of  the  disturbance  w.  Note  that  the  results 
in  Proposition  IV.2  is  not  applicable  to  the  system  in  (19) 
because  its  linearization  at  the  origin  is  not  asymptotically 
stable.  Nevertheless,  the  nonlinear  reachability  analysis  as 
outlined  in  Proposition  IV.  1  is  still  applicable. 


Proposition  IV.3.  Let  x\  £  7Zni ,  £2  €  Rf2 ,  and  w  £  Pnu’ 
and  consider 

ii(t)  =  Ari(f)  +  b(xi(t),x2(t))  +Ew{t) 
i2  (t)  =  q(xi  (t)) 

where  b  :  7^ni+Tl2  — ♦  1Zni  whose  entries  are  bilinear 
polynomials  in  xi  and  £2,  q  :  Pni  — ►  Rf2  whose  entries  are 
quadratic  polynomials  in  x\t  and  E  and  A  are  real  matrices 
of  appropriate  dimension  such  that  there  exist  Q\  =  Qj  >-  0 
and  e  >  0  with 

f  ATQi+QiA  QiE  1 

ETQi  -I  J  -  ~eL 

.  Then ,  there  exist  positive  definite  V  £  K[(£i  ,  £2)] ,  s  € 
E[£i],  and  R  >  0  such  that  bm  £  S[(£j ,  £2,  ^t;)]  where 

bm(x i,x2,w)  ■—  -  [Wf(xi,X2,w)  -wtw  +  {R2  -  V)s] 


Proof:  Let  V(x)  :=  xjQiXi  +  X2Q2X2.  where  Q2  = 
Q\  >~  0.  Then,  there  exist  B\,  B2,  Hi  >-  0,  and  H2  >~  0 
such  that 


xjQMxux*) 

xjQ2q(xi) 

xjQiXixJxi 

xjQ2X2X^Xi 


=  xf  Biz(xi,x2) 

=  xf  B2z(x1  ,  x2) 

=  z(xuxi)tMiz(xi,xi) 
=  z(xi ,  X2)TM2  z(ll ,  X2) 


and  —bm  can  be  decomposed  as 


£1 

T 

£l 

w 

D 

w 

Z(£j,£i) 

z(£i,£i) 

z(£j,£2) 

z(£j , £2 ) 

where  D  is 

A? Q\  +  Q\A  4-  ocR2 1  Q\E  0  Pi  4- P2 
ETQi  -10  0 

0  0  —a  Mi  0 

bJ^  ■+■  B2  0  0  — QA/2 

and  D  negative  semidefinite  by  proper  choice  of  a  (suf¬ 
ficiently  large)  and  R  (sufficiently  small).  Consequently, 
bm  €  E[(£1jX2,U;)].  □ 

V.  Region-of-attraction  analysis 

The  material  of  this  section  is  adapted  from  [24]  where 
similar  results  were  proven  in  the  context  of  robust  region-of- 
attraction  analysis  for  systems  with  parametric  uncertainty. 
For  simplicity,  we  focus  on  the  case  without  uncertainty. 
Consider  the  autonomous  nonlinear  dynamical  system 

*(*)  “  /(£(*)),  (21) 

where  x(t)  £  7Zn  is  the  state  vector  and  /  is  an  n-vector 
with  entries  in  R[x]  satisfying  /( 0)  =  0,  i.e.,  the  origin  is  an 
equilibrium  point  of  (21).  Let  <£(f;xo)  denote  the  solution 
to  (21)  at  time  t  with  the  initial  condition  £(0)  =  xo-  The 
region-of-attraction  of  the  origin  for  the  system  (21)  is 

(xo  €  7Zn  :  lim  <£(f;x0)  =  o\  . 

I  t— »oo  J 

A  modification  of  a  similar  result  in  [2]  provides  a  character¬ 
ization  of  invariant  subsets  of  the  ROA  in  terms  of  sublevel 
sets  of  appropriately  chosen  Lyapunov  functions. 

Lemma  V.l.  Let  7  £  R  be  positive .  If  there  exists  a 
continuously  differentiable  function  V  7Zn  — ►  1Z  such 
that 

fivvy  w  bounded ,  and  (22) 

V(0)  =  0  and  V(x)  >  0  for  all  x  £  lZn  (23) 
flvA  {0}  C  {£  €  Tin  :  W(x)/(£)  <  0}  ,  (24) 

then  is  an  invariant  subset  of  the  ROA.  <3 

In  order  to  enlarge  the  computed  invariant  subset  of 
the  ROA,  we  define  a  variable  sized  region  = 

{£  £  7Zn  :  p(x)  <  /?},  where  p  £  R[£]  is  a  fixed  positive 
definite  convex  polynomial,  and  maximize  /?  while  imposing 
the  constraint  QPyp  C  f along  with  the  constraints  (22)- 
(24). 

SOS  programming  and  simple  generalizations  of  the  S- 
procedure  (namely  Lemmas  ELI  and  II.2)  provide  algebraic 
sufficient  conditions  for  the  constraints  in  Lemma  V.l. 
Specifically,  let  l\  and  I2  be  a  positive  definite  polynomials. 
Then,  since  li  is  radially  unbounded,  the  constraint 

V  -lx  €E[£]  (25) 

and  V (0)  =  0  are  sufficient  conditions  for  (22)  and  (23).  By 
Lemma  II.l,  if  S\  £  E[x],  then 

-[(/?-  p)Si  +  (V  -  7)]  €£[*]  (26) 

implies  the  set  containment  fip $  C  and  by  Lemma 


II.2,  if  52,53  G  £[x],  then 

-  [(7  -  V)s2  +  W/s3  4-  h\  e  £[x]  (27) 

is  a  sufficient  condition  for  (24).  Consequently,  t  is 

a  subset  of  the  ROA  and  is  an  invariant  subset  of  the 

ROA,  where 

Proa  :=  max  0  subject  to  (25)  —  (27), 

ROA  V€V, v  1  v  '  (28) 

V(0)  =0 ,5i  G  £[x],/?>0 

and  V*  and  7*  are  optimal  values  of  V  and  7  in  (28).  Here, 
the  sets  V  and  Si  are  prescribed  finite-dimensional  subspaces 
of  polynomials. 

We  now  focus  on  systems  governed  by  ordinary  differen¬ 
tial  equations  of  the  form 

X  =  f(x)  =  Ax  +  f2{x)  +  /3(x), 

where  f2  and  /3  are  vectors  of  (purely)  quadr; 
cubic  polynomials,  respectively,  and  A  G  7£nxn,  and  prove 
that  asymptotic  stability  of  the  linearized  dynamics  is  also 
sufficient  for  the  feasibility  of  the  constraints  in  (28)  (for 
sufficiently  small  7  >  0). 

Proposition  V.l.  Let  f  be  an  n-vector  of  cubic  polynomials 
in  x  satisfying  /( 0)  =  0,  and  let  P  >-  0,  R\  >-  0,  i?2  0, 

p(x)  :=  xTPx,  /i(x)  :=  xTi?ix,  h(x)  :=  xTR2x. 

If  there  exists  Q  >-  0  such  that 

atq  +  qa*o, 

then  the  constraints  in  (28)  are  feasible  for  some  R  >  0.  < 

Proof:  The  proof  is  constructive.  Let  z  =  z(x)  be  as 
defined  in  Lemma  II.4,  Q  x  0  satisfy  ATQ  4-  QA  ^  — 2i?2 
and  Q  R\  (such  Q  can  be  obtained  by  properly  scaling 
Q).  Let 

e  :=  Amin(i?2),  V(x)  :=  xTQx , 

and  H  >-  0  be  such  that  (xTx)V'(x)  =  zTHz  (which  exists 
by  Lemma  II.4).  Let  M2  G  Rnxn*  and  symmetric  M3  G 
Rn,xn‘  satisfy 

VK/2(x)  =  XTM2Z 

W/3(x)  =  zrM3z. 


Define 


Sl(x) 

Cl 

=  v -f-jS 

_ 

—  Ar 

s2(x) 

=  C2XTX 

7 

_  € 

0 

=  _I_ 

23 1 

S3(x) 

=  1, 

(Mg'  +  yjM^  M2) 

• — wr — 


since  Q  -  Ri  ^  0. 


&i(x)  :=  -  [(7  -  ^52  +  VK/53  +  i2] 

1 T  r 

x 


where 
B\  := 

and 


— 7C2/  i?2  4-  QA)  —  M2/2 

-M2t/2  c2H-M3 


Bi  h 


§/  —  M2/2 

-Mf/2  c2H-M3 


>-  0 


by  the  Schur’s  complement  formula.  Consequently,  61  (x)  G 
E[x].  Finally, 


(30) 

□ 


(29) 

‘  1  ' 

T 

— 0s\  +7  0 

'  1  ‘ 

and 

X 

0  siP  —  Q 

X 

b2 


where  B2  y  0  and  consequently  62  G  E[x]. 


VI.  Interpretation  and  demonstration  of  results 

It  is  worth  re-stating  that  the  motivation  here  is  theoretical 
rather  than  practical.  The  conclusions  can  be  summarized 
as  that  the  nonlinear  local  analysis  (based  on  S -procedure 
and  SOS  programming  relaxations  as  proposed  in  [5],  [20], 
[7])  is  always  capable  of  returning  a  feasible  result  (i.e., 
corresponding  optimization  problems  are  feasible)  whenever 
corresponding  conditions  for  the  linearized  dynamics  are 
feasible.  Alternatively,  these  results  may  also  have  some 
limited  practical  value  in  constructing  (possibly  conservative) 
quantitative  results  for  the  nonlinear  system.  For  example. 
Propositions  V.l,  III.2,  and  IV.2  can  be  directly  used  to 
construct  feasible  solutions  for  the  problems  in  Eq.  (28) 
and  Propositions  III.  1  and  IV.  1,  respectively.  Proofs  of 
Propositions  V.l,  III.2,  and  IV.2  also  suggest  a  recipe  for 
constructing  less  conservative  feasible  solutions  for  these 
problems  by  searching  for  an  “optimal”  quadratic  Lyapunov 
function  (instead  of  fixing  V  to  a  Lyapunov  function  for 
the  linearization).  A  construction  in  the  case  of  region-of- 
attraction  analysis  can  be  summarized  as  follows:  Choose 
the  multipliers  5i,  52,  and  53  in  the  form  given  in  the 
proof  of  Proposition  V.l  with  the  free  parameter  C2.  Affinely 
parameterize  H,  M2>  and  M3  in  terms  of  Q  (note  that  there 
may  be  multiple  possible  parameterizations  for  M2  and  M3 
and  the  choice  may  change  the  quantitative  results  -  here 
we  arbitrarily  choose  one  parametrization).  Then,  nPyp*  is  a 
subset  of  the  ROA  where 


where  for  a  symmetric  matrix  M,  M+  denotes  its  projection 
into  the  positive  semidefinite  cone.  Clearly,  si  G  E[x],  52  € 
E[x],  and  53  G  £[x].  Note  that 

V(x)  -h(x)  =xt{Q-Ri)x  6  E[x], 


0*  :=  max  0  subject  to 

'r,C2,ff,Q=QT>:Ri 

—7 oil  -R2-  ATQ  -  QA  -Mi{Q)/2 

-M2(Q)t/ 2  ciH(Q)  -  M3(Q) 

-0  +  7  0  1 
0  P-Q  -U- 


h  0 


(31) 


Note  that  the  above  problem  can  be  solved  through  a  series 
of  convex  SDP  problems  by  a  line  search  on  C2.  Construction 
of  feasible  solutions  for  the  problems  in  Propositions  m.  1 
and  IV.  1  can  be  developed  in  a  similar  manner. 

The  value  of  such  “suboptimal”  construction  of  feasible 
solutions  for  the  problems  in  the  context  of  nonlinear  system 
analysis  may  be  better  appreciated  by  recalling  the  fact  that 
one  of  the  main  difficulties  in  SOS  programming  based 
nonlinear  system  analysis  is  the  computational  complexity 
of  the  SOS  programming.  The  procedure  outlined  above 
provides  an  ad  hoc  way  of  generating  (possibly  high  quality) 
solutions  for  the  corresponding  optimization  problems  or 
initial  seeds  for  further  optimization.  The  following  example 
demonstrates  this  construction  for  ROA  analysis  and  com¬ 
pares  the  results  with  “optimar’solutions  from  (28). 

Example  VI.l.  Consider  the  Van  der  Pol  dynamics 

±1  =  — X2 

i 2  =  H  4*  (i?  -  l)l2, 

which  have  a  stable  equilibrium  point  at  the  origin  and  an 
unstable  limit  cycle  around  the  origin  which  is  the  boundary 
of  the  ROA  of  the  equilibrium  point.  In  this  example ,  we  will 
construct  invariant  subsets  of  the  ROA  using  the  problem  in 
Eq.  (28)  and  Proposition  VI.  Let  x  =  [xj  ,  p(x)  =  xTx, 
li(x)  =  h(x)  —  10“6xTx.  The  solution  of  the  problem 
in  Eq.  (28)  with  a  quadratic  Lyapunov  function  candidate . 
(purely)  quadratic  S2,  and  scalar  s  \  and  S3  certifies  1.57 
to  be  a  subset  of  the  ROA.  The  feasible  solution  provided  in 
Proposition  V.  I  certifies  20  to  be  a  subset  of  the  ROA. 
Alternatively ,  by  the  procedure  outlined  above  certifies  that 
Op, 0.65  is  in  the  ROA.  < 

VII.  Conclusions 

Sum-of- squares  programming  based  analysis  of  nonlinear 
systems  with  polynomial  vector  fields  may  be  regarded  supe¬ 
rior  to  analysis  based  on  linearized  dynamics  in  the  sense  that 
the  former  is  capable  of  generating  quantitative  certificates 
as  opposed  to  conclusions  from  the  latter  valid  only  over  in¬ 
finitesimally  small  neighborhoods  of  the  equilibrium  points. 
However,  sum-of-squares  based  approach  involves  multiple 
relaxations.  Therefore,  it  is  not  obvious  if  the  sum-of-squares 
programming  based  nonlinear  analysis  can  return  feasible 
answers  whenever  linearization  based  analysis  does.  In  this 
paper,  we  proved  that,  for  a  restricted  but  practically  useful 
class  of  systems,  conditions  in  sum-of-squares  programming 
based  region-of- attraction,  reachability,  and  input-output  gain 
analyses  are  feasible  whenever  linearization  based  analysis  is 
conclusive.  Besides  the  theoretical  interest,  such  results  may 
lead  to  computationally  less  demanding,  potentially  more 
conservative  nonlinear  (compared  to  direct  use  of  sum-of- 
squares  programming  formulations)  analysis  tools. 
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